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Structures and Dynamics Laboratory 
NASA: Marshall Space Flight Center 
Huntsville, Alabama, 35812 
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Abstract: A review of the application of single particle hydrodynamics in models 
for the exchange of interphase momentum in continuum models of multiphase 
flow is presented. Considered are the equations of motion for a laminar, 
mechanical two phase flow. Inherent to this theory is a model for the interphase 
exchange of momentum due to drag between the dispersed particulate and conti- 
nuous fluid phases. In addition, applications of two phase flow theory to 
de-mixing flows require the modelling of interphase momentum exchange due to 
lift forces. The applications of single particle analysis in deriving models for 
drag and lift are examined. 
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I. Introduction 


Before the advent of numerical approximation techniques, coupled with 
large computational facilities the study of single particle hydrodynamics was 
popular. Following the early conceptual efforts on theories of multiphase flow 
[I], [2] investigations into the form of certain classes of phenomena were 
undertaken. Specifically, those investigating the form for models of interphase 
momentum exchange found a wealth of information by evoking the arguments 
once used in single particle hydrodynamics. As an example, consider the follow- 
ing quote [3], used in a context specific to arguments relating to forms of 
momentum exchange models in continuum theories of multiphase flows: 

"On the other hand, the terms ... and ... have no analogs in single particle 
calculations and will be neglected." 

It would appear that the application of single particle hydrodynamics in contin- 
uum models of multiphase flow has received a degree of acceptance. 


II. Continuum Theories of Multiphase Flow 

Continuum theories of multiphase flow have developed along parallel lines: 
Mixture Theory and the Theory of Interpenetrating Media with Moving Inter- 
faces. Both approaches are rooted in the same fundamental assumption, namely: 
That both the dispersed and continuous phases of the flow can be treated as and 
described within an Eulerian kinematic framework by the conservation equations 
of a macro-continua. Implicit in this assumption is that the variable fields of 
each phase are unique and continuous over the flow domain. The limits of this 
assumption for the case of dilute concentrations of the dispersed phase have been 
explored [4] and the alternative of a Lagrangian or ’particle tracking’ kinematic 
scheme for the dispersed phase forwarded. In addition, continuum models have 
been adapted to granular material flows [5] where the dispersed phase concentra- 
tion approaches a maximum. 

Mixture theories arise from the specialization of the classical field theory 
requirements of internally consistent thermodynamic arguments [6],[7]. In 
contrast, the theories of interpenetrating media address more directly modifica- 
tions to the classical transport equations due to discontinuous or ’jump’ conditions 
at moving phase boundaries [8],[9],[10],[1 1], To reduce to a local form, the 
conservation equations resulting from the theories of interpenetrating media must 
be averaged in either space or time. In fact, the key difference between the 
mixture theory formulation for multiphase flows and the averaged conservation 
equations for interpenetrating media is that; while the averaging process is 
implicit in the mixture theory approach it is an explicit operation in the course 
of writing the conservation equations from the interpenetrating, moving phase 
boundary approach. 

With few exceptions, it is reassuring to note both approaches result in the 
same set of conservation equations. 

Consider a simple two phase flow. That is one in which the dispersed 
phase is a dilute, mono-dispersed suspension of non-reacting, smooth, rigid, 
spherical particles in an incompressible, linearly viscous fluid. Both phases and 
the surroundings are thermally equilibrated. Laminar flow conditions prevail. 



Regardless of the method of formulation, the conservation equations reduce to the 
expressions for continuity and balance of momentum [10],[12] and constitute the 
continuum equations of motion of the mechanical theory of two phase flow. 

CONTINUITY: +5 • 4 , Vy" • « = D (DISPERSED), C (CONTINUOUS) ( 1 ) 

SATURATION: $°+<t» c =1 ^ 

BALANCE OF MOMENTUM: (^- + v“ y y“) = “ ( j ) n 5p n +5 l' I +«|> a p‘ , h"+A a + (x-p (3) 


WHERE $ (i , t ) : CONCENTRATION 
p : MATERIAL DENSITY 
y (X,t) : VELOCITY 

p (X.t) : SPHERICAL STRESS (PRESSURE) 

T(x,|): DEVIATORIC STRESS 

4 U.t) : INTERPHASE MOMENTUM EXCHANGE 
X (X.t) : SATURATION (CONTACT) PRESSURE 

The flow has two velocity fields and two concentration fields. Each phase has a 
unique, constant material density and may be acted upon by a set of external 
potentials or body forces. In addition, there exists unique expressions for the 
spherical and deviatoric elements of the dispersed and continuous phase stress 

tensors. Lastly, under the assumption that the flow is saturated, i.e. no phaseless 

voids may develop, the phases are coupled by an interface "saturation" pressure 
and momentum exchange or transfer between the phases. 

If the external potentials are specified, the momentum exchange terms 

modeled and the deviatoric elements of the stress tensors specified by constitutive 
assumption or neglected via arguments with respect to magnitude, the resulting 
is the unclosed system of 9 equations in 10 unknown fields; velocity, concentra- 
tion and spherical stress or pressure, for each phase. Commonly, the assumption 
of equal phase pressures is used in an effort to create a closed determinant 

system of fields and conservation equations [1],[11]. However, in should be 
noted that the wisdom of this assumption has been challenged on both physical 
and mathematical grounds [12],[13],[14]. 


ITT. Models of Momentum Exchange in Two Phase Flow 

As a common denominator, all theories of two phase flow embody some 
model for the exchange of momentum between phases. For simple, mechanical 
two phase flow it can be demonstrated rigorously that the sum of the interphase 
momentum transfer must be conservative, i.e. the sum of all momentum transfer 
terms is zero. Those exceptions to this "summing rule" are more a matter of 
bookkeeping than conceptual difference [15]. 



In addition, the degree of coupling between the phases of the flow both 
in a physical and mathematical sense is controlled by the momentum transfer 
model. Two way coupling, i.e. momentum transfer from one phase to the other 
and vice versa is implicit in the requirement of conservative transfer of momen- 
tum. In attempts to simplify the computational complexities associated with 
applying continuum theories of two phase flow the assumption of one way cou- 
pling is often evoked [16]. One way coupling allows for the transfer of 
momentum from the continuous to the dispersed phase, but not vice versa. In 
this case the process of momentum exchange is not conservative. Within this 
context, one way coupling is synonymous with the statement that the velocity 
field of the continuous phase is unchanged due to the presence of the dispersed 
phase within the flow. The computational simplifications result from the fact 
that it is no longer necessary to solve the conservation equations for the 
continuous phase field variables simultaneously with the conservation equations of 
the dispersed phase. Given a solution to the continuous phase variable fields, 
perhaps generated by single phase analysis or experimental techniques, the 
uncoupled conservation equations can be solved for the dispersed phase variable 
fields. 

The use of a step function to describe the "effectiveness" of momentum 
transfer has been proposed and applied to the problem of laminar two phase jet 
flows [1],[17]. This approach allows the degree of coupling to vary, step-wise 
from uncoupled to one way coupled. The arguments raised are: that in flows 
where interparticle spacing in large relative to the sum of the particle diameter 
and twice the fluid boundary layer thickness on the surface of the particles no 
net momentum transfer between the phases occurs. It is argued that when the 
interparticle spacings are large and the suspension is dilute, the slip between the 
phases results only in unrecoverable dissipation in the particle wakes. However, 
if the multitude of efforts in the analysis of two way coupled, dilute two phase 
flows can be used as an indication, then it would appear that the "ineffective- 
ness" of interphase momentum transfer in dilute two phase flows is not generally 
accepted. 

Arguments have been presented for the set of variables which constitute 
the general class of admissible momentum exchange functions [3],[18]. 


- <j> D s(v°-y c )+<t> D b 3-3 v° + c3(3-y c ) +l(v d v c Vd c + 

' ~ ( 4 ) 

WHERE: (*,,)= 1/2 (?Z Y<= + ^c j ) RATE AT DEFORMATION TENSOR 

where ... indicates that additional functions, of higher order in gradient, do exist. 
S,B,C and L are constructed from the scalar invariants of the admissible vector 
and tensor fields. In addition, owing to the discrete nature of the dispersed 
Phase, the admissibility, as a class of momentum exchange function of gradients 
of dispersed phase field variables is still debated [11]. If the existence of 
smooth, continuous first partials derivatives to the dispersed phase variables is at 
question, consider Implicit in the assumption of the continuum model of two 
phase flow was that the dispersed phase could be treated as a macro-continua, 
i.e. that the dispersed phase field variables are smooth and continuous! 
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In fact, the use of Divergence (Gauss’) theorem in the process of reducing the 
global conservation equations to the local form requires the existence of continu- 
ous first partial derivatives of dispersed phase velocity and concentration fields. 

The appropriateness of using phenomenologically based arguments for the 
purpose of identifying the specific forms of the momentum exchange models from 
the general admissible class is justified [19]. In fact, more often than not it is 
the lack of a phenomenological or physical analogue that results in the neglect of 
a class of the admissible momentum exchange functions. 

Depending on the purpose of the analysis, four generic categories of 
momentum exchange processes and models are identifiable: 


-Drag Forces 
-Lift Forces 

-Inertial Coupling or Virtual Mass Effects 
-Inertial History Effects 

The arguments for the inclusion of the inertial coupling/virtual mass 
effects and inertial history effects are raised during the construction of models 
for momentum exchange in two phase flows [ 1 1,(20],:2 1 ]. The analogous single 
particle hydrodynamic forces have been investigated and debated [22], [23], [24], [ :>]. 
Inertial coupling stems from the analysis of the forces required to 
displace a given volume of fluid during the acceleration of a particle through it. 
Likewise, the inertial history or Basset force is linked with the acceleration 
history of a particle moving through a quiescent fluid. However, the lack of 
complete agreement on the single particle hydrodynamic analysis of these effects 
and the lack of agreement on the forms or even the necessity for retaining these 
effects [12] in models of momentum exchange in two phase flows precludes them 
from additional discussion at this point. 


Drag Forces 

Common to all theories of two phase flow is a model for momentum 

exchange due to drag between the fluid phase and the particles of the dispersed 

phase [1 11,[12],[26],[27],[28]. Intrinsic to each of these models is the existence of a 
slip velocity between the phases. The result being a net drag force on each 

phase: 


AV'AV s0tD -^ 


where the factor of proportionality; S is derived from arguments which are 

rooted in single particle analysis. . 

The most sophisticated models for S are derived from an analysis of the 
mean, terminal sedimentation velocities of a dispersion of spheres falling through 
a quiescent Newtonian fluid under gravity [11],[12],[29]. The analysis is limited 
to flows, about any one single sphere in the Stokesian regime, where the Rey- 
nolds number between the fluid and the sphere is of order unity or smaller. 
This implies either very small slip velocities and/or a very viscous fluid. The 
net result being that the inertia of the fluid phase is neglected. 


Equilibrium between the force on a single sphere due to the gravitational poten- 
tial and the terminal or steady state sedimentation velocity drag can be given as: 


uo = ^(p°-p c )a < 6 > 

WHERE: a: PARTICLE RADIUS 

U: FLUID VISCOSITY 
a-' gravitational potential 

This sedimentation velocity potential is modified by the probability of hydrody- 
namic interaction with the next closest sphere, where the distance separating these 
spheres is defined probabilistically for a uniform dispersion. The result is the 
mean sedimentation velocity potential, corrected by the presence of a dispersion 
of other spheres of a given concentration: 


U - U 0 (1-6.55 $ D ) 


( 7 ) 


The incorporation of this result into a form for the interphase momentum 
exchange in two phase flow due to drag requires the suppositions that: At a given 
slip velocity the potentials acting on each phase due to drag are equal and 
opposite, i.e. momentum exchange is conserved. The slip velocity of two phase 
flow is then equated with the modified mean sedimentation velocity due to gra- 
vity. Resulting in, for a single spherical particle: 

- 6-55 ♦°) (v d -Y c ) - (p D - P c ) a 


This result is then further generalized with the assumption that if this is the 
potential for interphase momentum exchange due to drag on a single sphere, then 
the drag force per unit volume of mixed flow should be simply this potential 
normalized by the local volume fraction or concentration of the particulate phase: 


$Lg= -&d, g = + 6.55 <t>°) (V°-Y c ) (9) 


The caveats are obvious: The single particle analysis is grounded in the 
assumption of Stokes flow for the motion or slip of the dispersed phase relative 
to the continuous phase. Unlike the slip velocity of two phase flow, the sedi- 
mentation velocity at a given concentration is a constant, being the result of a 
constant potential, gravity. The modification to the mean sedimentation velocity 
by the hydrodynamic interactions with the other particles of the dispersion has 
presupposed that the dispersion is uniform, i.e. gradients of dispersed phase 
concentration are not present. 
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Lift Forces 


It is observed that certain classes of laminar shear flows will result in 
demixing of the dispersed phase [30], [31], i.e. the dispersed phase will become 
distributed in a non-uniform manner, despite the fact that initially the flow may 
have been homogeneous with respect to dispersed phase concentration. Analysis ot 
these phenomena have been made using single particle hydrodynamics [32] and 
using continuum two phase now theory with the 

description of the interphase momentum exchange [1 1],[12],[15],[1 oJ, 144J. it is 
generally accepted that the lift forces are those potentials, accounted for within 
the interphase momentum exchange which produce dispersed phase motions or 
migrations transverse to the slip or velocity difference between the phases. These 
lift forces arise through the interaction of the slip velocity, the rotation or spin 
of the particles of the dispersed phase and the shearing or gradients of the 
continuous phase velocity. If the exchange of momentum between the dispersed 
and continuous phase is conserved, then the lift force potentials act equally and 


in opposite sense on both phases. . 

A variety of single particle analysis, with the resulting identification o 

certain lift forces have been performed [23],[34],[35],[36],[37],[38]. The lift force 
on a single particle, as outlined by Saffman, [35] is the most commonly used m 
models of momentum exchange due to lift forces in continuum theories of 
demixing two phase flow. Saffman’s analysis identifies the lift force: 




( 10 ) 


WHERE v : FLUID KINEMATIC VISCOSITY 


u (y) : FLUID VELOCITY {IN PARALLEL FLOW) 

Unlike the assumption of Stokes flow used in the analysis of and subsequent 

application to momentum exchange due to drag, Saffman’s lift analysis begins 
with the Navier-Stokes equations. However, the Reynolds numbers for the sup 
velocity, the particle spin and the fluid shear are all constrained to order unity 
or smaller. Hence, even though the Saffman’s analysis retains the inertial terms 
of the Navier-Stokes equations, the flow is not inertially dominated. The 

solution requires the matching of the "inner" and "outer" asymptotic expansions 
of the flow equations in an analysis technique pioneered before numerical 
approximation coupled with computational methods became available [39]. In 

Saffman’s analysis the field variables for the flow about a single sphere are 

expanded about the radial position from the sphere center. The inner expansion 
has as a boundary condition the no slip requirement at the sphere surface Since 
the sphere is spinning, the no slip boundary condition constrains the fluid to 
have the angular velocity of the sphere surface. Hence, the particle spin enters 
the calculation implicitly, despite the fact that due to the level of truncation, it 
does not appear explicitly in the expression for lift. The necessity for an outer 
expansion results from the non-convergent nature of the solutions to the inner 
expansion as the distance from the sphere center approaches the infinite. 



The outer expansion embodies a second boundary condition, at a some large dis- 
tance from the sphere center; namely, the undisturbed (by particles) velocity field 
as radial position approaches the infinite. In fact the primary difference in most 
single particle analysis of lift is whether the boundary condition for the fluid 
velocity field in the outer expansion is constrained by a wall condition [32], a 
quiescent fluid [34] or by the rate of fluid strain [35]. The assumptions implicit 
in Saff man’s analysis of the lift force on a single sphere include: the flow is 
uniform and parallel, the slip velocity is parallel to the plane of the fluid shear, 
the shear or velocity gradients of the fluid are linear and the particle spin vector 
lies in the plane of the fluid shear, but is normal to the slip vector. The 
resulting force is normal to the plane of the fluid shear (and slip vector) as well 
as being normal to the spin vector of the particle. If, in terms of the slip 
velocity, the particle lags behind the fluid the lift will produce a migration of 
the particle into the faster, adjacent fluid and vice versa if the particle leads the 
fluid. In other words, the sense of the lift force depends on the sense of both 
the gradients of the fluid velocity and the slip velocity. 

Saff man s result for the lift on a single spherical particle has been 
generalized for the analysis of the momentum exchange due to lift forces in 
continuum theories of de-mixing two phase flow [1 1],[12],[15]. The lift force on 
the dispersed phase is normalized by the number of particles in a unit volume of 
two . phase flow. The resulting generalized form of the momentum exchange due 
to lift being written: 




<fr P c IP 


IY 


(ID 


This form of the momentum exchange due to lift has been used in calculations 
of parallel, de-mixing two phase flows using continuum theories. However, it is 
not clear that this general form will reduce directly to Saff man’s result for a 
one dimensional, parallel flow, both in terms of magnitude and the directional 
nature of the lift force. In addition, despite the fact that the intentions are well 
motivated, the meaning of operations such as absolute value, square root and 
division by a second order tensor valued variable is not clear. 


IV. Summary 

In conclusion, it has been shown that single particle hydrodynamics is the 
only source presently used to derive and justify forms for the interphase momen- 
tum exchange models within continuum theories of laminar multiphase. 

Four generic classes of momentum exchange models can be identified: 
Drag, Lift, Inertial or Virtual Mass effects and Inertial History or Basset forces. 
The latter two categories are still advent in nature and have not yet assumed a 
role in the models of interphase momentum exchange in applications of contin- 
uum theories of two phase flow. On the other hand, examples of Drag and Lift 
forces in applied momentum exchange models are numerous, though not without 
obvious caveats and inconsistencies. 

It is encouraging that the requirements of emerging technologies based on 
an understanding of multiphase flow processes has motivated such a great deal of 
work on generalized models for interphase momentum exchange. 
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Abstract 

A general framework is outlined for the modeling of fluid-particle flows. The 
momentum exchange between the constituents embodies both lift and drag forces, 
constitutive equations for which can be made explicit with reference to known single- 
particle analyses. Relevant results for lift are reviewed, and invariant representa- 
tions are posed. The fluid and particle velocities and the particle volume fraction are 
then decomposed into mean and fluctuating parts to characterize turbulent motions, 
and the equations of motion are averaged. In addition to the Reynolds stresses, fur- 
ther correlations between concentration and velocity fluctuations appear. These 
can be identified with turbulent transport processes such as “eddy diffusion” of 
the particles. When the drag force is dominant, the classical convection-dispersion 
model for turbulent transport of particles is recovered. When other interaction 
forces enter, particle segregation effects can arise. This is illustrated qualitatively 
by consideration of turbulent channel flow with lift effects included. 

Introduction 

Flow of an incompressible, single-phase fluid is fully characterized by a single 
kinematic field, the velocity. The kinematics of a fluid-particle mixture involves the 
velocity of each constituent, and an additional scalar field representing the particle 
volume fraction, or concentration. In some flows, the latter can be quite critical. 
For example, since the effective viscosity of a suspension is a strong function of the 
concentration, particle segregation in a viscometer will violate the assumption of 
homogeneity required to interpret measurements. Furthermore, if the distribution of 
particles is in part determined by the shear rate ( e.g ., Ho and Leal, 1974; McTigue, 
et a/., 1986), the apparent viscosity will be rate-dependent, and the mixture will 
appear to be non-Newtonian even when this may not be so locally. 

A great deal of work has been done on the dynamics of a single particle in a 
viscous fluid; reviews are given by Happel and Brenner (1965), Goldsmith and Ma- 
son (1967), Brenner (1966, 1970), and Leal (1980). In many applications, however, 
it is neither practical nor even of interest to track individual particles. Rather, the 

^his work was supported by Sandia National Laboratories under contract to the U. S. Department 
of Energy (DE-AC04-76DP00789). 
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primary concern is more often for some average characteristics of the flow. The 
objective of a continuum mixture theory is to provide governing equations for these 
average kinematic fields. Ideally, one would like to draw upon knowledge gained 
from single-particle analyses to guide the development of the constitutive models 
required by the continuum theory. As a practical matter, this process relies heavily 
upon empirical input as well. 

This paper is intended to illustrate by example the construction of a two-phase 
flow model for turbulent mixtures. It is highly idealized and far from complete, but 
captures some interesting phenomenology. We first outline the general mechanical 
balance laws for a mixture. We then review in detail results from the literature 
on lift forces in viscous and inviscid flow's. Generalizations in forms appropriate 
for the exchange of momentum between the constituents in a mixture are then 
discussed. “Exact” equations of motion are posed for the simplest forms for lift and 
drag interactions. Turbulent decomposition and averaging yields not only Reynolds 
stress terms, but other correlations of velocity and concentration fluctuations as 
well. Simple “eddy viscosity” and “eddy diffusivity” closure schemes are adopted 
to model the correlations. We then show that the classical convection-diffusion 
model for turbulent transport emerges naturally for the case when the drag term 
dominates the disperse phase momentum balance. Finally, we consider channel flow 
with the lift force present, and identify an equilibrium particle segregation due to 
a balance of lift and turbulent diffusion. 


Balance Laws 


The balance equations for the mass and momentum of constituent a are given 


by: 


dp a 

dt 


+ V ■ {p a V*) = 0, 


( 1 ) 


= V-T* + p a g + m 0 , (2) 

where p a is the density (mass of constituent a per unit volume of the mixture), 
v a is the velocity, T Q is the stress, g is the acceleration due to gravity, and m a 
is a body force due to the interaction of constituent a with the other constituents 
present. Equations (l) and (2) take the form of the classical balance laws for a 
single-phase continuum, with the exception of the interaction force, or momentum 
exchange, m a . Equation (l) neglects chemical interactions or phase changes, which 
would be embodied in mass exchange terms (c/, Passman, et a/., 1984). 

We anticipate that the mixture can be represented as a single continuum, so 
that: 

g + V ■ („v) = 0, (3) 


Pa 


' d\ c 
~di 


+ v a • Vv fl 


P 



+ v • Vv 


V • T + *>g, 


( 4 ) 
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where p, v, and T are the density, velocity, and stress for the mixture. Comparison 
of (l)-(4) shows that the mixture quantities are related to the constituent quantities 
by the summation rules : 

p = Zp. (») 

p\ = Ep a v a ( 6 ) 

T = E[T a -p a (v-v a )(v-v a )], ( 7 ) 

Em a = 0, (®) 

where E indicates the summation over all constituents present. Equation (8) shows 
that whatever momentum is lost from one constituent is gained by the other (s). 

The density fields, p Q , can vary due to changes in both the volume fraction, <j> a , 
and the local density (mass of constituent a per unit volume of that constituent), 
-y a . Thus, it is convenient to introduce the decomposition: 

Pa = 4>ala- 

Finally, we consider only saturated mixtures, in which all space is occupied, which 
imposes the requirement: 

E<£ a = 1. ( 1 °) 

Equations (l)-(lO) hold for multiphase systems with any number of constituents. 
For present purposes, we specialize to the case of two, a continuous fluid (a — J) 
and a dispersed particulate solid (a = s). In this case, we let 

d> = (f> s = l — <j>f. (n) 

We also restrict attention to mixtures comprised of incompressible constituents h/ 
and constant). 

Without loss of generality, it is convenient to decompose the stresses, T a , into 
an isotropic pressure, p a , and an extra stress , T a : 

T a = -(p a Pal + T*. ( 12 ) 

There is substantial motivation to include in the momentum exchange a buoyancy 
force, p/V<£, due to the fluid pressure acting over the interfacial surfaces (e.g.. 
Passman, et al., 1984). Thus, we also define an extra momentum exchange , m^, 
such that 

m, = p/V(^> + m* (13) 

Finally, equations (l)— (13) can be combined in the form: 

- ^ + v-K 1 -^) v /l = 0 ’ M 

ot 


d 4: + V • (*v.) = o, 


(15) 
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d\f 

dt 


+ Vf ■ Vv/ 


(16) 


P> 



v, • Vv, 


- -(1 - 4>)Vp f + V • T} + p f g - m*, 


~4>Vpj - V[<£(p, - Pf ) ] + V * T* 4 - P ? g 4 - m ;, 


(17) 


Note from (14) and (15) that, even though the constituents are taken to be incom- 
pressible, neither velocity field is, in general, divergence-free. 


Lift Forces 

Particle segregation has been observed experimentally in Poiseuille flow by a 
number of investigators; results have been summarized by Brenner (1966), Cox and 
Mason (1971), Goldsmith and Mason (1967), and Leal (1980). In general, these 
studies show that particles lagging the fluid motion tend to migrate toward the 
centerline, or region of minimum shear rate, and particles leading the fluid migrate 
toward the wall. Segre and Silberberg (1962) found that, for a small range of 
mean flow Reynolds number, neutrally buoyant particles can achieve an equilibrium 
position at a dimensionless radius of about 0.6. 

Saffman (1956) and Bretherton (1962) have shown that a particle embedded in 
a steady, rectilinear Stokes flow, i.e. at zero Reynolds number, cannot experience a 
net force normal to the unperturbed fluid streamlines. Thus, any analysis for the 
cross-stream lift on a particle in a steady, rectilinear flow must take inertia into 
account. One approach is to introduce small inertial effects through a perturbation 
of the Stokes flow problem, and a number of such analyses are in the literature. 
Both unbounded and bounded domains have been addressed. The analyses for the 
former assume that the boundaries of the unperturbed flow are sufficiently far away 
that they do not interact with the disturbance due to the particle. Boundaries play 
an indirect role in this type of problem, of course, insofar as their presence may 
be required to establish the velocity gradient or curvature with which the particle 
interacts. Two well-known analyses for unbounded flows are those by Rubinow and 
Keller (1961) and Saffman (1968), which are summarized briefly below. 

Analyses for bounded flows address configurations in which, say, a fixed wall 
lies within the disturbance field of the particle. Examples include the work of Ho 
and Leal (1974) and Vasseur and Cox (1976). It is not immediately apparent how 
one might adopt analyses of this type in the formulation of a continuum model. 
Although we have attempted previously to do so (McTigue, et a/., 1986) using 
the results of Ho and Leal, the result is not very satisfactory. The indication of 
difficulty is the appearance of the channel width in the expression for the lift. It 
would seem that this should enter through boundary conditions rather than through 
a constitutive equation. Obviously, this arises because the bounded-flow analyses 
are geometry-specific. For this reason, we consider in more detail generalizations 
only of lift forces in unbounded flows. 

Rubinow and Keller (1961) consider a sphere spinning with angular velocity fi 
and translating at velocity V' through an incompressible viscous fluid. The fluid is 


assumed to be static far from the sphere. The solution takes the form of a Stokes 
expansion in the near field that satisfies boundary conditions at the particle, but 
fails far away, and an Oseen expansion in the far field that exhibits the converse 
behavior. An asymptotic match is performed in order to calculate the forces on the 
sphere. The expansions are in powers of the particle Reynolds number, 


Ry 


7 jdV 


( 18 ) 


where a is the particle radius, V is the magnitude of the translation velocity, and 
H is the fluid viscosity. The result of interest here is that for the lift force normal 
to the direction of translation, ^ 


fj™> = tt aS/O x V'(l + O(Rv) 1- ( 19 ) 

Consider a rectilinear shearing flow, v fl {x 2 ). A force-free particle spins with the 
angular velocity of the fluid, so that H 3 = -k/2, where /c = dv fl /dx 2 is the shear 
rate, and V( = -V == -{v fl - v, 2 ). The “slip-spin” lift force (19) is then 

fg K > = \*o‘i,kV. ( 20 ) 

It is interesting to note that, although Rubinow and Keller’s analysis is for a small 
inertial correction to the Stokes flow problem, the lift force is, to leading order, 
independent of viscosity. 

Saffman (1968) considers a sphere in a simple shear flow, translating parallel to 
the undisturbed streamlines with a relative velocity of magnitude V'. The analysis 
again is based on matched asymptotic expansions. In addition to (18), two other 
Reynolds numbers enter the problem: 

= 2^, fl„ = 2^5, (2i) 

n M 

and the conditions under which the analysis holds are 

R v < R'J 2 , Rv < 1, An C 1- ( 22 ) 

For a simple shear flow given by v;\ = Vo + KX 2 i Saffman obtained a slip-shear 
lift force in the form: 

tjZ ) =6A6aW/V li [ sgn^Kr'V. (23) 

Equation (23) indicates that a particle lagging the fluid (F > 0) migrates toward 
higher-velocity streamlines, and a particle leading the fluid [V < 0) migrates in the 
direction of decreasing fluid velocity. This qualitative behavior is in accord^ with 
experimental observations. In the rectilinear shearing flow, the ratio of the “slip- 
spin” lift (20) found by Rubinow and Keller to the “slip-shear” lift (23) treated by 
Saffman, then, scales like R\'"- 


In a more general flow field, the applicability of the “slip-shear” analysis requires 
that the characteristic length scale for the disturbance field is much less than that for 
the variation in shear rate. Saffman (1968) suggested, from dimensional reasoning, 
that the lift due to interaction of the particle disturbance field with the mean-flow 
curvature takes the form 


f L 2 C) =: ca 4 ^ 2 f /3 n 1/3 K( sgn^)|ci 2/3 , (24) 


where ? is the curvature ( e.g ., for an undisturbed mean flow v fl = v 0 + kx 2 + 
fx 2 /2). Saffman noted that determination of both the sign and the magnitude of 
the constant c await more complete analysis. In this unbounded flow, one may 
define a curvature Reynolds number R ( = 7 /o 3 ^///, in terms of which the ratio of 
the “shear-curvature” lift (24) to the “slip-shear” lift (23) scales like R l J*R* ,3 Ry l . 
It is interesting to speculate upon the possibility that these two forces oppose one 
another in certain flows. Consider, for example, plane Poiseuille flow carrying a 
neutrally buoyant particle. The particle slips relative to the fluid due to the Faxen 
effect only. Thus, for symmetric flow in a channel of half-width d, the shear rate is 
K = ~Zvx 2 /d 2 , the curvature is 0 = -3 v/d 2 , and the slip is V = a 2 v/2d 2 , where v is 
the mean velocity. 1 he “slip-shear” and “shear-curvature” forces are then balanced 
where 


£2 

d 


0.80 c~ 2 R 1/3 , 


(25) 


where R = ^jvdjn is the channel Reynolds number. Segre and Silberberg (1962) 
observed an off-axis peak in particle concentration in flows of dilute suspensions in 
circular tubes. The peak occurred at a dimensionless radius of about 0.6, and was 
manifest in flows characterized by R of order 10. If these conditions apply to a plane 
geometry, the constant c would be of order 0.8. The ratio of the “slip-shear” (23) 
to the “shear-curvature” (24) lifts in Poiseuille flow scales like R~ 1/6 . Ho and Leal 
(1974) also considered interaction with the mean flow curvature, but in a bounded 
flow. The ratio of the curvature effect discussed by Saffman in an unbounded flow 
(24) to that found by Ho and Leal scales like R~ 1/3 . 

Both Rubinow and Keller and Saffman studied small Reynolds number effects. 
In the other limit, Drew and Lahey (1987) have recently considered mviscid rota- 
tional flow past a sphere. They obtain a lift force of exactly the same form as that 
found by Rubinow and Keller (19), but multiplied by a factor 4/3. The same result 
was obtained independently by Auton (1987). This is essentially like the classical 
Kutta-Joukowski lift on a two-dimensional body in a plane flow, which is just 7/LT, 
where U is the velocity of the body and T is the circulation. 


Invariant Forms for the Lift Force 

The momentum exchange, n^, includes fluid-particle interaction forces such as 
lift and drag. For brevity, let us decompose m; into drag, m* D , lift, m* L , and other 
components: 
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m* = m* D + m2 + . . . . 


(26) 



It has been suggested previously (Drew, 1976; McTigue, et al., 1986; Passman, 
1986) that the lift might include terms of the form 

= 2a 2 d>D / • (v y - v s ) + 4/? 2 d>D/ • (V • D/), (27) 


where D a = symVv a . It is expected that a 2 and /? 2 may be a functions of the 
particle volume fraction, <j>, the relative speed, |v/ - v.|, and the invariants of D/, 
D,, and their higher-order derivatives. In particular, if we assume that, for dilute 
suspensions, we should recover the single-particle results discussed in the foregoing 
section, 2 this function can be made explicit. For example, Saffman s result for the 
“slip-shear” lift (23) is recovered for the choice 


3(6.46) / 7/M 2 \ 1/4 
a2__ 4 jra \2trD 2 / 


(28) 


The “shear-curvature” lift (24) is recovered for the choice 


_ 3 ccnf/ /* 1/3 (29) 

132 ~ 47r|2V-D / i 1 / 3 ' 

Generalization of a “slip-spin” lift of the form found by Rubinow and Keller 
(20) poses some difficulty. It would appear that such a lift is proportional to 
2W; • (v/ - v,), where W, = skwVv/ is the skew-symmetric part of the fluid 
velocity gradient. However, W/ is not invariant ( e.g ., Truesdell, 1977, p. 115). 
Drew and Lahey (1987) have suggested that this dilemma can be resolved by si- 
multaneous consideration of the virtual mass effect. The virtual mass, too, when 
generalized from the classical expression, is not easily put into an invariant form. 
However, the combination of the virtual mass and lift forces posed by Drew and 
Lahey is invariant: 


* 



P/vA 
Dt ) 


— 2W / • (v/ — v s ) , 


(30) 


where my M ’ s t ^ ie momentum exchange due to the virtual mass effect, and the 
substantial derivative is defined by D a /Dt = d/dt + v a ■ V. In (30), neither the 
virtual mass, represented by the difference in convective accelerations, nor the lift, 
in the form 2W, • (v/ - v,), is invariant, while their sum is. This depends upon 
the remarkable result that the coefficient lf /2 is the same for both the virtual mass 
and the lift. That (30) embodies the classical virtual mass effect is easily seen by 
specializing to an unsteady, uniform flow. That it embodies the result of Drew and 
Lahey for the lift can be demonstrated by specializing to steady, rectilinear shearing 
flow. Drew and Lahey point out that a simple regrouping of terms can yield an 


invariant form for the virtual mass: 

* 1 A( D>yf} 

mvM = fit* H-dT " 



V,) • V(v/ - v.) , 


(31) 


and a lift in the form of the first term in (27) with a 2 - O'// 2 - Equations (31) and 
the lift sum to recover (30). 


2 This assumption was stated by Drew (1976) as the principle of correct low concentration limits. 



Turbulent Decomposition and Averaging 

It is evident from the foregoing discussion concerning lift forces that the for- 
mulation of the necessary constitutive models necessary to complete the equations 
of motion (14-17) is quite formidable. It remains to specify relationships for the 
stresses, T} and T*, the pressure difference, p s - p f , and momentum exchanges 
such as that due to drag, m* D . Each of these raises subtle and complex modeling is- 
sues. For present purposes, we skirt these difficulties in order to isolate phenomena 
associated with the lift and drag. In particular, we assume 


T / = (32) 

T ** = 0, (33) 

Pf=P, = P, (34) 

m z> = «i<£(v/ - v,), (35) 

= 2a 2 <t>Df • (v, - v,). (36) 


We rationalize neglect of the fluid extra stress, T} (32), by confining attention to 
inertially-dominated flows. 3 In a dilute suspension, it is easy to imagine that the 
disperse-phase extra stress, T a *, vanishes (33), implying that there is no direct ex- 
change of momentum between particles. The assumption of equal pressures (34) 
implies that Brownian motion (Nunziato, 1983) and certain inertial effects at the 
particle scale (Givler, 1987) are negligible. The drag force, m^, is written in its 
familiar form (35), proportional to the relative velocity. The coefficient Oj is, in 
general, expected to depend upon 4> and [v/ - v,|, accounting for the effects of 
particle interference at high concentration and inertia at high relative velocity, re- 
spectively. The choice «i = 9p/2 a 2 corresponds to the classical result for Stokes 
drag on a single particle, and is adopted here. Finally, we take the lift in the form of 
(36), with q 2 assumed to be constant for simplicity. Neglect of inertial effects in the 
drag (35) while retaining those giving rise to the lift (36) is justified if RvR~ n « 1 , 
where n = 1 for the “slip-spin” lift of Rubinow and Keller (20) and n - 1/2* for the 
“slip-shear” lift of Saffman (23). 

Under these assumptions, the momentum equations (16-17) reduce to: 


p f {jot + V/ ' Vv/ ) = -(l-^)Vp + p/g 

-a.i<t>{\ / - v a ) - 2a. 2 <t>T)f • (v/ - v.„), 
+ V * ' = -‘fr'Vp+Psg 


+ a i H y f - v ») + 2«2<£D/ • (v/ - v s ), 


(37) 


(38) 


3 T^e dissipative, viscous terms represented by T} are critical, of course, to the extension of this 
discussion to the kinetic energy balances. 



Note that these “exact” equations of motion do not contain interaction terms rep- 
resenting diffusive forces. 

The independent field variables in (37) and (38) are each decomposed according 
to: 


a — v a + v l, 

(39) 

4> - <(> + 4>\ 

(40) 

p = p + p, 

(41) 


where overbars indicate mean quantities and primes indicate fluctuating quantities. 
By definition, — 0 and p 7 = 0 . However, the averaging scheme chosen here defines 
the mean velocities in terms of mean momenta , an approach introduced originally 
for compressible, single-phase flows (Favre, 1965), and suggested in the multiphase 
context by Drew (1975): 

P a v a = P^Va- ( 42 ) 

Note that p a = 7 a <f> a for incompressible constituents. Note also that the averages of 
the velocity fluctuations do not vanish, but (1 — </>)v’y = 0 and — 0. 

Substitution of (39)-(41) into (14), (15), (37), and (38) and averaging yields: 


d4> 

dt 


+ V • [(1 - *)v/] = 0 , 


d<}> 

dt 


+ V • (tj>v s ) = 0 , 


(43) 


(44) 


P f f ^ + v/ • Vv/j = -(1 - $) Vp + + V ■ T'y + p f g 

-<*1 [0(vy - V,) + 4>v' f ] 

-2 a 2 \d> D/ • (V/ - v,) + D/ ■ <^Vy + ■ (v, - v,)] 

(45) 

p, (^ + v,- Vv,) = + + 

+ [0(v/ - Vj) + <p\' f ] 

+2a 2 U D/ • (v, - v a ) + D/ • W f + ■ (v/ - v,)] 

(46) 

where = “PaV^v' a is a Reynolds stress for constituent a, and triple correlations 
have been omitted. 



Constitutive Models for Turbulent Correlations 

The averaged mass balances (43-44) appear in forms identical to the exact equa- 
tions (14-15), in part as a consequence of the definition of average velocity (42). 
Averaging the momentum balances, however, yields a number of correlations of 
fluctuating quantities. The Reynolds stress terms are familiar from single-phase 
turbulence, but several additional correlations arise here that are a direct conse- 
quence of the fluctuations in the particle concentration field, <t>. Of special note is 
the correlation of particle concentration and fluid velocity fluctuations, which 
represents a flux of particles due to the fluid turbulence. In all subsequent devel- 
opments, we neglect the pressure-concentration correlations appearing in (45) and 
(46). 

For the present discussion, it suffices to adopt the simplest possible closure 
scheme, following essentially the classical “eddy viscosity” argument. That is, a 
correlation of some fluctuating quantity with v' a is Jaken to be proportional to the 
gradient of the mean of that quantity: 


- PaV' a v' a - u a l lQ [V ■ (p a v a )]l + 2u a / 2Q sym V (p Q v a ) , (47) 

<t>\' } = -u f l 3 V$, (48) 

2 = -u / / 4 V(V^), (49) 

where u a is an appropriate velocity scale for constituent ct. and the Is are appropriate 
length scales (“mixing lengths”). 

Turbulent Convection and Dispersion of Particles 

A commonly encountered situation for which modeling capabilities are well de- 
veloped is that for particles fully entrained in the fluid. In this case, the particles 
are essentially passive’ tracers for the fluid, and are transported by the mean con- 
vective motion and by turbulent diffusion. It is worth considering briefly where this 
classical model is embedded in the mixture theory outlined here. 

For <f> 1, the fluid mass and momentum balances (43 and 45) are approxi- 

mately those for the fluid alone: 


Vv 7 = 0, (50) 

If + = -Vp + V • + 7/g. (51) 

In this approximation, the fluid motion is unaffected by the presence of the particles, 
and can be solved independently. The disperse-phase mass balance (44) remains in 
its exact form. Suppose the drag coefficient, Qj, is large, so that the dominant terms 
in (46) are simply those due to drag. This can always be realized for sufficiently 



small particles; the ratio of the drag to lift forces discussed in the foregoing scales 
at least like a” 1 . The disperse-phase momentum balance then reduces to. 

4>{\j - v ? ) = ( 52 ) 

t.e., the mean flux of particles relative to the fluid is balanced by the turbulent 
correlation </>vy. 

Equations (44), (50), and (52) combine to give. 

d A + y f . V^= -V • {W})- ( 53 ) 

dt 

Substitution of (48) into (53) yields: 

d A +V/ . = V-(PV^), ( 54 ) 

dt 

where V = u// 3 . This recovers the classical result: the particle concentration field is 
governed by a convection-dispersion equation, with turbulent dispersion coefficient 
or “eddy diffusivity” D. A similar discussion for the case when the gravitational 
body force is retained was presented by McTigue (1981, 1983). 

Although this limiting case is relatively simple and quite well known, the present 
development is revealing. Many texts derive the turbulent diffusion equation solely 
from a statement of mass balance, a decomposition and averaging process, and a 
model for the correlation of concentration and velocity fluctuations. This tends to 
mask the fact that the turbulent diffusion is a dynamic process in response to fluid- 
particle interactions. Thus, the momentum equations must be considered. Indeed, 
it is worth reiterating here that the turbulent dispersive flux, appearing in in 
(53) arose from decomposing and averaging the drag force (35). Thus, the tendency 
for the particles to be convected with the mean fluid velocity and dispersed by the 
fluid velocity fluctuations is clearly identified with the drag. Analogous observations 
have been made previously with regard to molecular diffusion {e.g., .Muller, 1968). 

We also note that the assumptions leading to (54) are quite special, and empha- 
size in particular the neglect of any interaction forces other than drag in writing 
(52). It is evident that much more complex phenomenology could be embedded in 
the general scheme outlined here if additional interaction forces come into play. 

Channel Flow With Lift Effects 

Particle segregation has been observed in turbulent jets (Laats and Frishman, 
1970), and ascribed to the “Magnus” lift force (30). Here we consider plane channel 
flow in order to simplify the kinematics, and retain the cross-stream lift effects 
embodied in (46). The analysis is highly simplified and somewhat speculative, and 
is intended only to illustrate the type of phenomena that might be represented by 
a mixture model of the type sketched out here. 



Consider a vertical channel, with steady upward flow. The flow is in the 
direction, so that g = {- g , 0, 0}, x 2 = 0 at the wall, and x 2 = h at the centerline. 
We assume that 7, > 7/, so that gravitational settling will cause the disperse 
particulate phase to lag the fluid. 


For a steady, rectilinear flow in the mean, the mass balances (43, 44) are identi- 
cally satisfied. We expect again that, for ^«1, the fluid momentum balance can 
be approximated by that for the fluid alone (51). Thus, in the streamwise direction, 
(51) becomes: 


dp 


d 




(55) 


For a smooth-walled channel, familiar arguments for the mixing length l 2f (47) and 
the identity u f = u* = (r 0 /7/) 1/2 , where r 0 is the shear stress at the wall, lead to 
the usual logarithmic velocity profile: 


V J1 

u* 


1 . U+X 2 

= - In 


k i / 


+ 5.5, 


(56) 


where k — 0.4 is the Karman constant, and v — is the kinematic viscosity. The 
streamwise momentum balance for the particles, from (46) and (48), and neglecting 
T * a 21 , becomes: 


0 ~ ~ ls< ^ g + ~ a 2U+K s x 2 ^-~, (57) 

ax i ax 2 ax 2 

where we have assumed Ufl 3 — u*/c 3 x 2 . Substituting from (55) for the mean pressure 
gradient in (57), and noting that the fluid shear stress gradient is simply — lf u l/h* 
(57) becomes: 

° = ~ 1 ^ 9 + ~ U *i) - • (58) 

Let us suppose, again for simplicity, that the second and third terms in (58) domi- 
nate. In this case, we are left with a balance between the buoyant weight and the 
drag, giving 

Vfi ~ v*i = Vco, (59) 

where Voo = 2a 2 (7^ — ')f)g/9/j. is the Stokes settling velocity. 

The cross-stream momentum balance for the disperse phase is 

/ 1 d<f> \ dv f 1 

0 = - 5 .,), ( 60 ) 

which is simply a balance between turbulent diffusion and the lift due to the mean 
flow. Note from ( 60 ) that the gradient of 4> vanishes at the centerline if the fluid 
velocity gradient vanishes there. According to the eddy viscosity model adopted for 
the Reynolds stress (47), the latter is in fact required by symmetry. However, of 


course, the logarithmic velocity profile (56) does not satisfythis condition. There- 
fore, we can anticipate a similar failing in the solution for <j>. Substitution of (56) 
and (59) into (60) and integration gives 


<j> = 4>(h) exp 


OClVocU* / 

l-r? y 

olikDk V 

V ). 


(61) 


where r? = x 2 /h, and a diffusivity, V h = u+K s h has been introduced. This profile 
has some of the expected characteristics: the particles are concentrated toward the 
center of the channel by the lift; the central peak is flattened by diffusion; and the 
channel margins, where the fluid velocity gradient is steepest, can be essentially 
clear of particles. That (61) indicates tf>(0) = 0 is a result of using the logarithmic 
fluid velocity profile (56), which is not valid in the limit i 2 -» 0, in (60). 

Lee and Durst (1982) conducted experiments in this configuration using glass 
beads in an air stream. Some of their results are in qualitative agreement with 
those found here: the air velocity profile (56) is little affected by the presence of 
the particles (at less than 0.5% mean volume fraction); the particles lag the fluid 
approximately by their fall velocity (59); and there is a particle-free zone near 
the wall. However, important phenomena are missed by this simple analysis. In 
particular, Lee and Durst observed that the velocity difference (59) is not uniform 
across the channel, but typically decreases toward zero near the wall. The particle 
velocity profiles, then, are more nearly uniform across the channel, suggesting that 
the turbulent mixing brings high-momentum particles from the core of the flow 
toward the boundary. For the smaller particles examined, the profiles actually cross 
near the wall; i.e., the particles lead the fluid, so that the momentum exchange due 
to drag (35) changes sign. These phenomena are clearly not embodied in the model 
analysis outlined here. The limitation is most likely in the simple, Boussinesq closure 
scheme adopted (47-49). Kashiwa (1987) has modeled these experiments using a 
higher-order ( k - c) closure, and is able to represent the cross-stream transport of 
streamwise particle momentum into the near-wall region. 


Summary and Discussion 

Treatment of turbulent suspensions in the context of the continuum theory of 
mixtures is currently in its most rudimentary stages. The appeal of the overall 
approach is that it provides an axiomatic framework on which to build. In practice, 
of course, one is quickly confronted with the difficulty of posing specific constitutive 
equations for the stresses and momentum exchange that embody the phenomena 
of interest. This is only compounded in the case of turbulent mixtures, in which 
correlations between the three kinematic fields, V/, v s , and (j> proliferate. The intent 
of this paper is not to lay out a definitive set of equations of motion for such a system. 
Rather, we have attempted only to sketch the general spirit of the approach, and to 
illustrate by means of the simplest possible example. The sequence is familiar from 
its antecedents in classical, single-fluid flow: state balance laws, pose constitutive 
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equations, construct exact equations of motion, introduce a decomposition and 
averaging scheme, model the resulting correlations, and, finally, solve boundary 
value problems. The channel flow problem considered here, exhibiting particle 
segregation effects, only hints at the rich and complex phenomenology that could 
be embedded in such a model. 

Each section of the paper encounters challenges. Exact forms for lift forces, 
even from single-particle analyses, are not well established; those that are known 
are complex, and their generalizations are not immediately obvious. We emphasize 
in particular that bounded flows have been analyzed ( e.g ., Vasseur and Cox, 1976) 
in which wall effects are critical, and it is not clear how one might adopt such results 
in a continuum model. Many of these remarks carry over to other interaction forces, 
as well, such as the Basset” term (e.g., Hinze, 1975, p. 463), which accounts for 
the history of the particle acceleration. Because no universally valid expressions 
for lift, drag, or other forces are available, considerable judgement is required in 
selecting the forms appropriate to a particular application. Constitutive equations 
for any concentration beyond the dilute limit are especially difficult to define; few 
analytical results are available (e.g., Batchelor’s (1972) work on “hindered settling”) 
and resort is usually made to empiricism. 

Perhaps the greatest challenge encountered in constructing a model for turbu- 
lent mixtures is the “closure” problem, familiar from single-phase turbulence, but 
magnified here by the presence of additional fluctuating fields. As in single-phase 
problems, some simple configurations can be addressed through classical Boussi- 
nesq models (e.g., 47 49) and simple scaling arguments. However, it is also clear, 
even from the highly idealized channel flow problem addressed here, that such an 
approach is severely limited. For example, the Boussinesq model for the Reynolds 
stresses (47) does not embody normal stress effects in rectilinear flows, while one 
might easily imagine that such effects could be important. Higher-order closure 
schemes are obviously called for, and steps in this direction have been taken with 
sorne success. Scheiwiller (1986) has developed a h — c model to represent snow 
avalanches, and has achieved excellent agreement with laboratory experiments. 
Kashiwa (1987) has used a similar approach, and successfully captures some of 
the unusual phenomena observed by Lee and Durst (1982) in the vertical channel 
flow discussed in the foregoing section. 
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INIHDCUCT1CN 

Gas -sol id two -phase flows occur commonly in many natural and industrial 
situations. Examples are blood flows, rocket exhaust plumes, pulverized 
coal gasification and combustion, and sediment transport by air and water. 
These flows are invariably turbulent and are characterized by the mutual 
coupling between the solid particles and the gas phase. Contrary to passive 
additives in a single -phase flow, the particles will change the flow 
structure of the carrying fluid. Globally, metering and heat transfer data 
[1,2,3] of two-phase flows shows discrepancy from the single-phase data. 
Further, small scale turbulence structures are also affected. Solid 
particles may attenuate the spreading rate and damp the turbulance intensity 
in a jet flow [5,6] . The alternation of the turbulence structure was found to 
depend on the particle size, the solid loading ratio as well as the physical 
properties of the different existing phases. 

In general, a complete theoretical treatment of two-phase flows is not 
possible because of the lack of detailed understanding of the physical 
processes involved [7] . Previous analytical studies have not been very 
successful, due in part to a lack of knowledge about the turbulent flow field 
of the conveying gas which is a prerequisite to the solution of the two -phase 
flow problem. Difficulties in theoretical analysis also arise from the 
coupling between the two phases, i.e., the exchange of momentum, mass and 
energy between phases. These coupling phenomena comprise a very complex 
interaction which affects both the gas and particulate phases 
Consideration of the infinite variety of interfacial geometries and flow 
regimes, various forms of non-equilibrium, and aggregation of particles 
complicates the problem even further. 

The inability of the theoretical analysis to account for all the 
complicated interactions in two-phase flows is similar in the study of 
single-phase turbulent flows two decades or so ago. An exact theory of 
turbulence did not (and still does not) exist; however, using a combination of 
theoretical equations, modeling assumptions, and experimental evidence, 
mathematical models describing certain features of the flow were developed. 
The field of turbulence modeling has subsequently been developed to the point 
where single -phase turbulent flow fields can be predicted rather well using a 
variety of turbulence models of varying complexity [8,9]. These advances 
suggest that a similar combination of theory, experiment, and modeling could 
be used to develop computational models capable of predicting two -phase 
flows. However, extra sets of equations and correlations need to be 
formulated and modeled for turbulent gas -solid flows. The purpose of this 
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paper is to discuss and review the recent advances in two -phase turbulence 
modeling techniques and their applications in various gas -solid suspension 
flow situations. In addition to the turbulence closures, heat transfer 
effect, particle dispersion and wall effects are partially covered here. 


APH5QAOES AD GCHERNDC HJKTKKS 


Because of the intrinsic, complex coupling between different species in 
two-phase flows, there seems to be no "unified" set of governing equations 
that can completely describe the flow field of two-phase media. However, 
there are quite a number of different formulations in the literature from 
which to begin. One approach, the so called "discrete" or "tracking" 
approach, starts with an equation of motion for a single discrete particle in 
a turbulent fluid flow field and the particle's trajectory is calculated. 
For particles much smaller than the smallest scales (say Kolmogorov's 
microscale) of turbulent motion and for which the solid's material density is 
much greater than the conveying gas, the BBO equation (Basset, Boussinesq, 
Oseen) , which is the momentum equation of a single particle, can be reduced to 
[ 10 ]: 


d 1 

— Vi = — (Uj - v i ) + gj 

dt t* 


( 1 ) 


Because the Eulerian velocity u i is a stochastic quantity when the conveying 
gas flow is turbulent, this simple looking ODE cannot be solved analytically 
due to its inherent nonlinearity. 

However, progress has been made using this approach in conjunction with 
the turbulence closure models which have been developed for single phase 
flows. The basic strategy is to use the turbulence model to calculate the 
fluid flow field assuming that no particles are present. This calculation is 
used to generate the velocity in equation (1) after making suitable 
assumptions regarding turbulent time scales, length scales and isotropy. To 
account for the mutual coupling (or the "two-way" coupling [11]) of mass, 
momentum, and energy between phases, the extra source terms generated by 
particles must be included in the Eulerian sets of governing equations for the 
gas phase. In the mean flow fields, this can be achieved by the iterative 
PSIC (particle-source-in-cell) technique developed by Crowe and co-workers 
[12,13] or by the non-iterative, transient numerical scheme of Dukowicz [14] . 
The discrete particle approach can also be extended to account for the 
particle- turbulence interactions which have two aspects -- the turbulent 
particle dispersion (the influence of fluid turbulence on the particles) , and 
the "modulation" effect [15] (the effect of particles on fluid flow 
turbulence). These will be discussed in further detail in the next section. 

In the non- discrete (continuum) approaches, two formulations are 
commonly used; the first considers the gas -solid suspension to be represented 
as a single inhomogeneous medium. The interactive forces between the phases 
are taken account of by internal stresses which must be related by 
constitutive equations to the bulk properties of the medium. Sets of 
governing equations for this approach were first formulated by 
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Barenblett [16] and described in detail in Monin and Yaglom [17]. This 
approach was also used recently in heat transfer analysis of a gas -particle 
pipe flow [18]. 

The other approach is the so-called "two -fluid” approach. This 
approach regards the gas and particles as two inter-penetrating continua in 
much the same way as the two species of a flowing binary mixture, for example. 
Here, the cloud of particles is regarded as a continuum and the governing 
equations are obtained by properly averaging the conservation equations over 
a volume and expressing the equations in differential forms. Many authors, 
namely Murray [19] , Drew [20] , Marble [21] and Ahmad i [22] , have described the 
two-phase flow based on the two -fluid formulation and applied it to some 
physical processes. It is often not possible to formulate a general set of 
governing equations for gas -solid two -phase turbulent flows due to the lack 
of understanding and differences in interpretation of the physical processes 
involved (for example, the "solid-phase pressure” term [23]). In order to 
obtain theoretical relations of two-phase turbulent flows, several 
assumptions have to be invoked to simplify the formulation. These are: 

1. The particle phase is dilute (volume fraction of particles, 

0 « 1) and is made up of particles spherical in shape and uniform in 
size. The oarticle material density P 9 >> P } so that the model is valid 
when p p = 0 ( p) This assumption is required because we 

ignore particle-particle collisions, the frequency of which increase 
quadrat ically with loading. The uniformity of particle size reduces 
the book-keeping in the formulation; extension to poly-dispersed non- 
uniform size distribution is a straight forward matter for dilute 
suspensions . 

2. Both the particulate and fluid phases behave macroscopically as 
continua. The fluid phase is Newtonian and both phases have constant 
physical properties and do not undergo any phase change. The continuum 
hypothesis assumes that the mathematical "points" are large enough to 
contain many particles and fluid molecules to ensure a stationary 
average. In order to satisfy the "dilute suspension" and continuum 
assumptions simultaneously for particle phases, some stringent 
restrictions regarding the number of particles in a smallest control 
volume made up of Kolmogorov microscale , the distance between particles 
to avoid direct inter-particle interactions, have been discussed 
[10,24,25]. However, the continuum approach has proven to be 
applicable also to situations which do not strictly meet such conditions 
[26]. 

3. The mean flow is steady and incompressible. Molecular diffusion, 
Brownian motion and gravity effects on the particulate phase are 
negligible compared with turbulent diffusion. Electrical and magnetic 
forces are not considered here. 

With the above assumptions, we may adopt the governing equations 
developed by Marble [21,27] and Hinze [25] which are applicable to dilute gas- 
particle flows. Marble used statistical averages for the particle cloud and 
postulated the macroscopic governing equation for the gas phase. Continuity 
equations are written for each phase: 

dp 3 

— + ( pU i ) = 0 ( 2 ) 

at dx { 
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9p p d 

+ (p pVj) = 0 (3) 

at axj 

Here P p is the mass of particles per -ini t volume of mixture (or "density" of 
the particulate phase where Pp = P e & / P is the particulate phase volume 
fraction). The momentum conservation equations for each phase are: 

3uj duj 3p 8 

P + pUj = ( r j j ) + F pi (4) 

at axj axj axj 

avj avj 

Pp + PpVj — = - F P i (5) 

at ax j 

Here F p - x represents the force acting on the primary fluid per unit volume due 
to the presence of the particle. Note that due to dilute assumption, the 
multiplication of (1 - <f> ) by each term inequation (2) and (4) was replaced by 
1. Of special note is that the values of P» P that appear in the 

continuum relations (2) and (4) are, in a sense, "smoothed” variables. The 
detailed gas disturbance caused by the particle motions are omitted from the 
instantaneous gas velocity vector Uj. Since the gas velocity varies 
strongly in the neighborhood of a particle that is moving through the gas, use 
of these smoothed variables in continuity, momentum and energy relations 
requires that all particle wakes or regions of immediate influence are 
dissipated very rapidly over the gas control volume. Hinze [25] treats this 
problem by attributing the forces around the particle as the external forces 
and disregarding the modified velocity field around the particles. If this 
external drag force follows Stokes law, then the fluid velocity in the 
Stokes drag law is at "infinity”, i.e., a large distance from the particle 
center so that the detailed fluid motion in the neighborhood of the particle 
is still not accounted for. However, the inadequacy of this model is not 
important for small volume fractions of particles having a not too large 
velocity relative to the gas. But for large volume fractions and cases in 
which particles may form into groups by trailing another in its wake, large 
regions of the flow may be inadequately modeled (c.f. [25] and [28]). 

Several derivations concerning the two -fluid model equations have 
appeared in the literature. The derivations include those of Hinze [29] , Soo 
[30], Drew and Segel [31], Ishii [32], Nunziato [33] and more recently, Roco 
and Shook [34]. The resulting equations differ in various ways such as the 
pressure gradient term for both gas and particulate phase, momentum source 
term, and shear stress tensor of the secondary phase although general 
constitutive equations relating stress and flow properties have not yet been 
developed. However, for the low concentration limit of suspension flows of 
small spherical particles, most of the derivations will recover similar 
forms. For example, the theory proposed by Ahmadi [35] which was general to 
the extent that it could be applied to both concentrated and dilute two-phase 
flows could be shown to recover the theory of dusty gas as derived by Saffman 
[36] in the low solid volume fraction range. The general expression for the 
internal forces between solid and continuous phase is discussed by Truesdell 
and Toupin [37] (also see [38]) and Drew and Segel [31]. The philosophical 
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reasons for using the two- fluid continuum approach and the common feature of 
dispersed two -phase flow systems can be seen in Drew's [20] review paper. 

By performing the Reynolds decomposition and time averaging of 
equations (2) - (5), the following mean equations for statistically steady 
flows result. 
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(7) 


( 8 ) 


F P l 


( 9 ) 


where some use of the continuity equations has been made in deriving 
equations (8) and (9). As a consequence of dilute suspension, triple 
correlations involving fluctuations in the particulate phase density are 
considered negligible. At this point the mean interaction term F . needs to 
be specified. Empirical expressions for the interaction terms nave been 
summarized by [39] for low and moderate solids concentrations. The 
appropriate relationships are given by 


__ 1 

F = — (l + 0. 179 -J Re p + 0 . 013Re p ) p p (V j - Uj) (10) 

t* 

+ _ (1 + ( 3 / 2 ) 0 . 179^| Re p + 0.113 x 2Re p ) p p * (Vj • - Uj 7 ) 

” here | U , - V j | d, 

Re p = (11) 


and t* “ <Jpp,/18^ (12) 

We note that several turbulent correlations appe ar _i n equations (7) 

(9) in addition to the conventional Reynolds stress u i v i for the single- 
phase flows. These terms arise from the velocity fluctuations and 
fluctuality volume fractions of the particulate phase and represent the 
turbulent momentum flux and mass flux of particles. To close the set of 
governing mean equations, models are required for these second order 
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correlations. The field of turbulence modeling for single-phase flows is a 

rapidly expanding one and will form the two-phase closure models described 
here . 


THD-H1ASE TIFBULENCE HEELING 

hierarchy of turbulence closure models has been received by Reynolds 
[40], and recently in [8,9], The proposals range in complexity from zero 
equation models where the turbulent fluxes are modeled as if they were 
molecular fluxes, with an eddy diffusivity related to mean flow structures to 
Mean Reynolds Stress models where separate transport equations are solved for 
each component of the turbulent flux vectors and tensors. Most two-phase 
turbulence models follow the single-phase turbulence models for 
incompressible flows closely; their modeling is discussed in the following. 

The most common and simplest modeling technique is to assume a Newtonian 
type constitutive equation for relating the turbulent fluxes to the mean 
field through an eddy viscosity. For gas phase Reynolds stresses, 


U i 'Vj ’ = - K fSij + / 3 a ijk 

1 

where S j j is the mean rate of strain tensor of gas flows — 


k = »/ 


2 U i u i and ? f is the eddy viscosity. 


(13) 

auj auj 1 

3xj axj ) 


Depending on the level of complexity employed, the eddy viscosity could be 
specif ied by zero -equation mixing length models , one-equation models or two - 
equation models. Due to the presence of solid particles, the eddy viscosity 
constructed by these models must take into account this effect. 


ZERO -EQUATION MODELS 


Early theoretical studies [41 - 43] indicate that the presence of solid 
particles decreases the eddy viscosity of the gas flows arising from 
dissipation of turbulence energy at the interface between solid particles and 
the fluid. These results lead several first-order closure schemes which 
modify the eddy viscosity for the clean gas flow without suspension of solid 
P ar ticles, v f 0 . For example, Owen [41] proposed 


( 14 ) 
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for t e /t* > 1 

This model has been used by Melville and Bray [44] for application in a 
turbulent free jet of dilute gas -particle mixture and has been further 
modified by Choi and Chung [45] and Chung et. al [46] for application in a 
wall-bounded shear flow. Most closure models developed at this level 
heavily involve empirical information and limiting case (loading ratio 
approaching zero) analysis and, in most cases, ignore the effect of particle 
size [47,48]. This level of models are very useful in engineering analysis 
because of their simple forms. However, they fail to handle some important 
effects, such as the "turbulence modulation", and besides, it is hard to 
prescribe the "mixing length" and the effect of particle on the mixing length 
scale. The next level of models, which incorporate a transport equation for 
the turbulence kinetic energy, and thus the velocity scale, were developed in 
the hope of providing additional generality and at the same time account for 
the effect of particles on the turbulence structure. 


ONE- EQUATION MODEL 

An equation describing the dynamics of the gas-phase turbulence kinetic 
energy can be derived from equations (4) and (8) by simple manipulations . 
For statistically steady, high Reynolds number flows it is given 

8k 9Ji 

Uj +P k - « + u i' F P i (16) 

8x j 9x j 


Here 
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is the rate of production of turbulence energy and 
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is the rate of energy dissipation rate. 
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and 


Uj ’Fpj the extra particle production (or dissipation) term, all per unit 
of mass. Probably the first attempt to use a one-equation turbulence model 
to study the two-phase flows is that of Dannon et. al [49], They applied a k-i 
closure model to a particle- laden axi- symmetric jet. The length scale SL was 
specified algebraically and was taken to be the same as that of a single-phase 
jet. For the k-equation (16), the diffusion term and production term were 
modeled following the conventional single -phase gradient- type modeling 
technique. In their study, quasi -equilibrium (i.e. Uj t Vj ) and 

mono-dispersed particles in Stokes regimes were assumed, 
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which simplify the interaction terms F p j 


Pp(Vi - Uj) 


and the turbulence 


kinetic energy equation becomes 
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The triple correlation was neglected and the concern was the modeling of the 
additional dissipation term created by the particle slip velocity at the 
fluctuation level. This term is similar to the turbulence "modulation" 
effect attributed to the inability of dispersed-phase particles to 
completely follow turbulent eddy fluctuations at high frequency. This added 
dissipation mechanism has been experimentally observed [2,5,6,50,51] and has 
gained much attention in recent two-phase modeling studies. 

The fluctuating velocity correlation of this term is bounded by 


0 - “ (U i 'Vj 1 - Uj ‘\l i •) < 2k ( 18) 

where the two bounds represent the cases where particles completely follow 
the fluid ( u ', = v () and stationary particles relative to the velocity 
fluctuation ( V £ = 0). Dannon et. al [49] proposed a model that has the 
correct limiting behavior 


- (Uj'Vj* - u j ' u i 1 ) = 2k [1 - exp (-B(t*/T))] (19) 


1/2 

where! = (We) is the 


Kolmogorov time scale, 


and 6 is a model constant. 


The use of the time scale T was argued [52] to be inappropriate since the 
eddies contributing most to the correlation are the energetic eddies 

which have an integral time scale t e . Dannon et. al [49] indicated that this 
model did not give good prediction due to a change in the structure of the 
turbulence and the structure was represented by the length scale. As a 
result, they had to arbitrarily modify the production and dissipation terms 
reflect the structural variations. Due to the difficulty of specifying 
the length- scale distribution a priori in a flow and appropriate modeling for 
particle effect, most workers have abandoned one -equation models in favor of 
two-equation or even stress-equation models in which the length scale is 
computed from a transport equation. 


40 



TWO- EQUATION MODELS 

Most studies in two-phase turbulence modeling utilizing a transport 
equation for the turbulence length scale if are based on a modeled equation 
of the isotropic dissipation rate £ ; this equation can also be derived from 
equations (2), (4), (8) by appropriate differentiation, multiplication and 

averaging. The exact £ -equation consists of 67 terms with particle's effect 
accounted for [53] . For high Reynolds number flows and based on an order- of - 
magnitude analysis [54] , the groups representing the production of £ by 
vortex stretching, the viscous destruction of £ , and the diffusive flux of £ 
in the * i direction, which are not affected by particles , are usually modeled 
following the single-phase k— 6 model of Jones and Launder [55] . The 
resulting modeled equation (except the extra particle destruction term) 
becomes 
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The last term in the RHS of the above equation is the contribution from inter - 
phase transport and is another main effort in two-phase modeling. This 
equation is solved simultaneous with the k -equation to estimate the eddy 
viscosity Y t = C u k 2 /*. Since the effects of particles are modeled through 
the k and * equations, is assigned 0.09 same as the single phase flows. 

In the two -equation model level, there have been several proposals for 
modeling the extra terms in the k and e equation. Chen and Wood [52] 
basically followed [49] and proposed exponential forms for added dissipation 
terms in both the and equation. In the k equation, the correlation 

U V | * Is modeled as 

uY T V r T r = 2k eX P (~ B k (21) 

where t Is the time scale of the energetic eddies and in the context of the 
model is given by kand € . here has often been interpreted as the 

lifetime of a typical turbulence eddy by [56,57] in their Lagrangian 
calculations. The time ratio t*/t € is the Stokes number [12] measuring the 
response of how quick the particle responds to a typical eddy turnover. To 
generalize the model, the constant Bfc was Introduced and was determined by 
limiting behavior for small particles, which corresponded to the linear 
perturbation analysis with respect to a passive additive by [42] . A similar 
approach has been used by Pourahmadi and Humphrey [58] and Gavin et. al [59] in 
the k- equation. Their model for this term is summarized: 

U^vP = 2k / (1 + t*/t e ) (22) 

It can easily be shown, for small values of that this model yields the 

same results as equation (21), depending on the numerical value of B^andthe 
form assumed for t e . Genchev and Karpuzov [60] assumed t* >> t e so that 
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u i’ v i* s Oand no model was required. However they also assume that the 
particles will follow the mean fluid motion which implies << t <t This 
inconsistency plus the lack of model comparisons with data in their paper 
casts doubt on their model. 

The model of Elgobashi and Abou-Arab [53] for the correlation u ( 'v t • 
based on Chao's [61] solution of the linearized Lagranglan equation of motion 
of a spherical particle in a turbulent flow. Their model, in its most general 
form, is extremely difficult to implement, especially for wall -bounded two- 
phase flows, owing to the necessity of computing definite integrals over all 
possible frequencies of fluid motion at every grid point. Given the 
uncertainties of the model it seems more appropriate to use a relatively 
simple model which exhibits the correct asymptotic behavior such as [52] and 
[58] (c.f. [62]). 

As in modeling single-phase flows it is the e -equation which provides 
the greatest uncertainty. In the model of Chen and Wood [26], particle- 
hydrodynamic drag force was assumed to follow Stokes law, thus the last term 
of equation (20) became 

T ? i • 3F. j • 2 ji \ aiii • ( avj' auj • \ ' 

2 ¥ v - (23) 

p 3x k 3x k t* p [ dXj l dXj dXj 


A similar exponential model was proposed for this term as given by equation 
(19). A different time scale, i.e. Kolmogov's time scale T, was used here in 
place of t # since the eddies contributing most to the high frequency 
destruction mechanisms are the dissipative eddies which have a time scale 


<»/«> 


9uj • #v t • 

T 3 ■ In [26] it was assumed that and are completely 

»x, ax , 


uncoupled on this time scale since >> ^ in most practical eas -solid 
turbulent flows. In this limit equation (23) becomes (2(/t*) (p^/ p) 
Clearly this model is not correct for very small particles or if z & ( 
In the model of [58] an additional term is added but it assumes that the added 


sink of dissipation to be a function of the integral time scale t # and hence 
does not seem to be particularly appropriate. The extra sink and destruction 
of ( are modeled collectively by [53] as C,j * t (*/t)in equation (20) where ( t 
equals ( plus the extra dissipation terms appearing in the k- equation and 
C*i is kept the same constant as in the single-phase flows. 


IHSEHtSHN CF SD5ESDBD SOLID BNCEKLES 

Extremely small particles which behave like trace molecules can be 
treated as "passive" contaminants in the turbulent flow field. The behavior 
of clouds of particles may be extended from single -particle dynamics when the 
mixture is very dilute, say, the volume fraction of solid < <9 ( 10" 2 )(c . f . 
[7]). The subject of passive additive transport has been treated 
extensively in the text of Monin and Yaglom [17] . See also the book by Hinze 
[63] and the review paper by Launder [64]. However, as the particle size 
increases, dispersion will be opposed by particle inertia and so once some 
critical particle size is exceeded, discrete particle dispersion must be 
treated in a different way from "passive" contaminant diffusion. 
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For dilute suspensions, the particle trajectories can be calculated y 
the tracking approach. In this approach, particle dispersion ue to 
turbulence has been modeled by random walk [65] or a Monte-Carlo Stochastic 
method [14 66]. These methods usually require extensive computational 
storage and time to achieve a stationary average In some Lagrangian 
approaches, certain types of diffusional velocity have been modeled for the 
particle motion which is usually proportional to the concentration gradient 
[67 68]. In the two-fluid approach, particle dispersion due to turbulence 
L represented by the correlations and/or . J* lo ™f„ le 7 el 

closures, these are usually modeled as a gradient type, Fickian diffusion 

process : 

- — = D 111 (24) 
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This constitutive equation is arbitrary at this point, however it may be 
justified theoretically under certain conditions [17] . 

Values of D t and D p are calculated by the value eddy v ^ COS J, tV _ / 
most models by introducing che turbulent Schmidt numbers^ Thus * 
and D, - r,/Sc, 0 r ’ V 8c » . where 'p is the effective eddy 

viscosity of the dispersed phase (to be discussed latex) This tjrpeof 

phenomenological approach for diffusion process However 

classical theory of "fluid point" diffusion of Taylor (c f . [63] ) . However 
when particle size increases, discrete particle diffusion is °PP osa y 

particle inertia and the crossing-trajectories effect [6 9 .7°]. Sin 

heavier particles have the tendency to "fall out" from one eddy to another 
the correlations between particles and fluid velocities decrease The 
effect is to decrease the particle dispersion. The effect of the pa 
inertia is not that clear. The inertia effect is characterized by the 
particle relaxation time t* , which is controlled by the physical properties 
of particles and the fluid, and the flow characteristics There have b een 
arguments concerning the characteristic flow time scales [71,72 73] f 
turbulent dispersion, although it has been indicated that the 
the heavy particles is a little larger than that of the light particle 
Higher-order modelings such as the one developed by [74] are able to predict 
this behavior . In the context of the lower -order phenomenological technique 
just mentioned, the Schmidt number should be modeled taking into accoun 
these effects. Recently one such model has been developed including a 
constant drift velocity [75] 


D 0 = 1 / (1 + 0.3 


|U, - V, 


/ V, ' v j ') 


1 / 2 


(25) 


in which the coefficient 0.3 was tuned based on the lateral dispersion of 

solid particles and measurements of [71,73]. 

However, in most models [26,53,58] the turbulent Schmidt number was 
simply set to some constant value following the turbulent mass transfer of a 
passive additive [76]. Some experimental data for gas-solid jets [77] 
however indicate that a constant value of Sc t (i.e. independent o g) 

is appropriate. For axi-symmetric flows, Sc t Sc p 0.7 is used [5 , 

In [53 58] the turbulent Schmidt number was simply cnosen to be one. 


43 



Finally, due to the continuum formulation, the correlation v i ’ v i ' has to 
be modeled. Following the gradient type model for the gas phase, this 

turbulent stress in the particulate phase is modeled from the Boussinesq 
assumption: M 
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Earlier analytical studies [79,80] have contributed to the understanding of 
some of the basic mechanisms of indirect interaction between particles 
through the surrounding particles. This has lead Melville and Bray [44] in 
their zero-equation modeling to propose 



(i + 


(27) 


Similar models have been used by [46,52] although the evaluation of t is 
somewhat different. Alonso [81] reviewed some developments in determining y 
and recommended the use of Peskins [82] formula p 


P 2 

— = 1 - (T L e/15 y) 
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whenK 2t*/T L and T L ; k/e. This model has been used bv [ 58,83], 
Although this model recovers the correct form in the limit t -»o (i.e. r - r 
as equation (27) , it will yield negative value of v p for reasonable values ol t * 
and T L taken from pipe flow data. It is not likely that** < 0 is 
physically appropriate, casting doubt on this model. 

The modeling technique discussed above based on the continuum approach, 
especially the modification of k and e equation, has been extended and 
adopted in some Lagrangian formulations. It has been shown by Shuen et. al 
[6] that using the stochastic formulation instantaneous properties are 
known; therefore, the extra dissipation term due to particles in the k 
equation is exact and requires no modeling. This calculation is rather 
complex and recently Mostafa and Mongia [75,91] have utilized the continuum 
two-phase model of [58] to model this term in Lagrangian calculations. A 
similar approach has been taken by [84] which highly simplified the 
stochastic calculations. The extra sink term in the (-equation is not 
closed in Monte-Carlo stochastic formulation and is modeled through a 
gradient type model in [6], The sensitivity of this term has been tested 
recently [85] and it has been found that this term is important in conjunction 
with the modulation term in the k-equation. Incorporation of the effects of 
particle on turbulence scale and the response of the dispersed phase in two- 
phase turbulence models is essential for representing the structure of 
particle -laden turbulent flows. 



tftLL EFFFCES 


For wall-bounded flows, boundary conditions for both gas and 
particulate phases are required. This is particularly important when other 
transport processes such as heat transfer and erosion are involved. The 
effect of particles on boundary layer and viscous sublayer flows has been 
studied analytically [10,86] . In most calculations, it is assumed that the 
influence of the particulate phase on the velocity defect law is to modify the 
logarithmic law in sublayer. Based on the Monin-Obakhov similarity analysis 
for the analogous stable stratified atmospheric boundary layer a set of 
"wall functions" taking into account the effects of particle size and loa ing 
was used by Chen [87]. The logarithm profile was modified as 
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with 6 ~ 5 and particulate flux Richardson number 
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The wall shear stress of the two-phase flow is related to that of the single- 
phase flow by (c.f . [10]). 


= 1 + 
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The wall boundary conditions for particulate phase are complicated by 
the unsteady particle-wall interactions such as deposition of particles on 
the wall and re -entrainment mechanism. The resulting piece of information 
from a multitude of particle -wall interactions can be assessed as the slip 
velocity for the particulate phase at the wall [88,89]. An expression for 
the slip velocity and wall shear stress was suggested by Soo [88] based on 
rarefied gas -dynamics theories. The lack of particle -particle interaction 
in accordance with the dilute suspension assumption gives a wall slip 
velocity : 
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and stress 
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where the fluid-particle interaction length £p*is given 
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3 , 1/2 

l P [ • “ Vj * ) 2 J t* 

and V y is set to be zero at the wall following an impermeable wall condition 
tor particles. These expressions have been utilized by [58 87] for their 
wall-bounded calculations. The effect of the walls on the particle drag 
coefficient, the particle turbulent intensity and the correlation Ui 'Vi ' 
and Ui * (U| * - V| • ) , which represent the particle-gas interaction and 

contribute to the modulation of the wall turbulence structures has been 
?^ died r ® centl y b y Risk and Elghobashi [90] by including Magnus’ force and 
lifting force in the particle dynamic equation. Their analysis can be 
incorporated for detailed gas -particle wall function development. 


SQMHDT AID nrsnKST rw 


Recent developments of two -phase turbulence models were reviewed. 
Most existing models are constructed following the familiar form of single- 
phase turbulence models. As in single-phase problems, most models are 
addressed through classical Boussinesq gradient-type diffusion processes 
and scaling arguments . Most models are also developed based on the treatment 
of turbulent suspensions in the context of the continuum, two-fluid theorv of 
mixtures . J 


The appeal of the two -phase closure technique embedded in the two -fluid 
continuum formulation is that it provides an axiomatic approach on which the 
analogous single-phase turbulence models are built. In practice, one is 
confronted with the difficulty of constructing specific constitutive models 
tor the stresses and momentum transfer, turbulent mass fluxes and mass 
transport in which additional fluctuating fields are magnified by the 
presence of solid particles. Following the modeling approach in single- 
phase flows, for simple flows (such as free shear flows) most two-phase models 
were addressed through classical Boussinesq assumptions and characteristic 
scaling arguments. Depending on the relaxation time scales, particles not 
only influence the higher wave number end of the gas -phase turbulence 
spectrum [10,25], but also the energy -containing range of the turbulence 
spectrum which is largely responsible for mixing [72,73]. Most proposals 
for treating turbulence modulations based on the two-equation k- ( model were 
not particularly successful for complex flows since they did not incorporate 
the turbulence scale effects and the response of the dispersed phase. 
Higher-order closure schemes or closures involving multiple- scale 

characterization of the gas turbulent spectrum are obviously called for, and 
some steps in this direction have been taken recently [74,87] Additional 
measurements similar to that of [72,73] are also needed to gain a better 

understanding of particle- turbulence scale interactions and modulations in 
multiphase flows. 

The real challenge and difficulty in developing a two-phase closure 
model for particle dispersion in connection with the two-fluid formulation is 
encountered in wall-bounded flows and poly-dispersed systems. 

Establishment of wall boundary conditions for the particle concentration and 
velocities depend on the interaction of particles with the wall. Particle- 
wall collisions are not always elastic and the phenomena is unsteady. The 
slip velocity boundary condition based on rarefied gas dynamics concepts is 
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probably not appropriate since the normal component of averaged velocity is 
not zero. Detailed measurements [92,93] and analysis [93,94] are needed. 

To extend the two- fluid continuum mixture theory for poly-dispersed 
situations, the continuous droplet models such as the one described in [95] 
can be used. The particles are represented by a statistical distribution 
function in a multi -dimensional space of droplet size, velocity location and 
time. The properties of particles are determined by solving the 

conservation of the distribution function. 

One merit of the two -phase turbulence models developed on the continuum 
formulation is that they can be accommodated into the Lagrangian approach and 
do not require excessive computational storage and time [75,84], 
Incorporation of more physics as turbulent combustion, evaporating sprays, 
boundary layer dust ingestion, poses no conceptual difficulties. Further 
testing and validation through well-defined experiments for more complex 
flows to establish the universality of the model constants are highly 
recommended. 
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ABSTRACT 


A multiple-realization particle trajectory scheme has been developed and applied 
the numerical prediction of confined turbulent fluid-particle £>«. The exmnpte 
flows investigated include the vertical pipe upflow experimental data of Tsuji et al. 
an°d trexperimenta, data of Leavitt for a coaxial jet now. oomp 
ticle-laden central jet and a clean annular je , » 

chamber. The results obtained from the numerical scheme agree 
experimental data lending confidence to the modeling approach. The mult ‘P le 
realization particle trajectory turbulent flow modeling scheme 15 belie ^ 
be a more elegant and accurate approach to the extension of single-particle 
hydrodynamics to dilute multi-particle systems than the more commonly employed 
two-nuid modeling approach. It is also better able to -corporate additional 
force terms such as lift, virtual mass and Bassett hrstory terms drrectly -to 
the oarticle equation of motion as appropriate. This makes . 

didate fo! particle migration studies and an extension to situations rnvo Ivrng 
liquid particulate phases with possible propulsion appl.cat.ons, such 
spray combustion, follows naturally. 
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INTRODUCTION 


anniirl» l - rbulent / 1Uld P3rtlCle f,ows are encountered in numerous technological 
comh T° nS SUCh „ aS fluidize d-bed combustors and pulverized coal gasifiers and 
combustors as well as in atmospheric studies involving the dispersion of pollutants 
The modeling of such turbulent flows involving the presence of a dispersed Z e 
made up of small, light particles further complicates the already complex phenomena 
encountered in single phase turbulent flows. However, the need to optimize t^TsLn 
process in technological applications involving turbulent fluid-particle flows or 

" a a " Ce d * . e pred,ct,on accuracy of atmospheric dispersion models makes it impossible 
to avoid the quest for a deeper understanding of the fundamental problems BesWel 

of fin 10U ff interaCtin8 C ° mplex P henomen a encountered in the modeling of this class 
of flows offer a very rich source of challenges to the fluid flow researcher 

i n The propu ! s,on systems for space transportation vehicles, in particular the 

L:r n 7 u fr L ' ieW; '"- 

na.n 8 ^ """ “ 

aung the elfects of the breakup or coalescence of droplets and bubbles from 

other particle-turbulence interactions encountered in such flows. 

two 7.° common a PP r oaches adopted in the literature for the modeling of 
two-phase flows are the homogeneous and the separated models. The former is 
applicable to situations in which the mean slip between the phases is small and the 

mirflu P x a es a Tn e s rS t f ? ° f ^ bUlk V3riety SUCh 35 the Pressure drop or 

n r S Wh6re m ° re det3i,ed inf ormation about intra- or inter-phase 

a ior is of interest, or where there is substantial segregation of the phases the 

iTz 

the Eulerian approach or whether a Lagrangian descriotion of th» a- S " ,b * d us,ng 
be more appropriate. 8 descr, P t,on of the dispersed phase will 

In the following, we present a discussion of turbulent fluid-particle flow 
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modeling in which the continuous phase is described using the continuum Eulerian 
approach while a Lagrangian description is adopted for the dispersed phase. We s a 
restrict ourselves to confined flows and thus include a discussion of the treatment of 
solid boundaries using the Eulerian - Lagrangian scheme. 


2. PARTICLE TRAJECT ORY SCHEMES 


In the Eulerian - Lagrangian modeling of two-phase flows, the continuous fluid 
phase is described using the standard single phase continuum equations. How- 
ever, the dispersed phase is modeled by computing for individual particles the 
trajectories and temperature histories where appropriate. The dispersed phase 
velocity and temperature fields are subsequently obtained from information 
obtained from the realization of a sufficiently large ensemble of particle tra- 
jectories. 

The use of a particle trajectory scheme in the modeling of turbulent fluid- 
particle flows represents only a subset in the field of computer simulation using 
particles as discussed by Hockney and Eastwood [1981]. Other important applications 
particle schemes discussed by Hockney and Eastwood include the modeling of covalent 
and ionic liquids, stellar and galaxy clusters, plasma and semiconductor devices. 

In fluid dynamic applications, the Particle-In-Cell (PIC) method of Harlow 
[1964] and, later, the Particle-Source-In Cell (PSI-Cell) method of Crowe et al. 
[1977] have received considerable attention. In the present investigation, the PSI-Cell 
method has been adopted as the basis for the Eulerian - Lagrangian model developed. 

The usual starting point for the development of fluid-particle flow theory is the 
consideration of the motion of a single particle in an infinite fluid. The nature of 
such a single-particle flow has been investigated by numerous researchers including 
Bassett [1888], Boussinesq [1903], Oseen [1927], Tchen [1947], Corrsin and Lumley 
[1956] ,Hjelmfelt and Mockros [1966] and Maxey and Riley [1983] and is relatively 
well understood for flows both within and outside of the Stokes flow regime. In the 
Eulerian treatment of the dispersed phase, the single particle flow theory is adopted 
directly to describe a multi- particle system and the validity of such a step is assumed. 
However, in the Lagrangian particle tracking approach, the focus remains on single 
particle hydrodynamics for obtaining an ensemble of statistical realizations, in this 
case the particle trajectories, which are then analyzed using the well established 
mathematical theory of statistics to extract the required phase information. 

In the presence of turbulence, particle trajectories are not deterministic due to 
an imposition over the mean velocity of a rapidly fluctuating random velocity 
component. This additional velocity component due to turbulence enhances the 
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dispersion of the particles, in aggregate, while the presence of the particles in the 
continuous phase, even in relatively small concentrations [Al-Taweel and Landau 
1977], does modify the underlying turbulence appreciably. This ’two-way cou- 
pling’ between the turbulence and the particulate phase exercises considerable 
influence over the evolution of such flows. These important effects will be 
considered later. 


3. GOVERNING EQUATIONS 


The field equations for the continuous phase in the Eulerian - Lagrangian 
scheme are the same as those for single phase flows except for the addition of an 
extra source term which accounts for the influence of the particulate phase on the 
continuous phase. The equations are written in a generalized form as 



+ S + Sp 


(3.1) 


where U] are the instantaneous velocity components, F e ff the effective exchange 
coefficients, S the usual single- phase source terms, S p the source terms due to the 
particulate phase and <|> any of the field variables such as velocity component, 
temperature for flows involving energy exchange, turbulence kinetic energy or its 
dissipation rate. 

The simplified form of the particle trajectory equation in which only the 
hydrodynamic drag term between the phases is retained [Adeniji-Fashola and Chen 


dvj (Uj-Vj) 

dt = ~ (3.2) 


where, m general, the fluid and particle velocities, uj and Vj respectively are made 
up of a mean and a fluctuating component and Z* is a particle response time 
defined in terms of the particle relaxation time t* which is valid for particle motion 
within the Stokes regime. Thus we have 


Uj = Uj +Uj' 


(3.3) 


v i = 


Vj +v/ 


(3.4) 
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X. = t./f 


(3.5) 


where 


t.= 


p sdp 
18 n 


(3.6) 


and 


CpRep 

24 


1 +0.15 Re r 


0.687 


FOR Rep = 


lUj-Vjldp ri 


FOR Rep > 1 


(3.7) 


An expression similar to equation (3.2) for the particle temperature history 
written for a particle thermal equilibration time tjjj can a i so be written for 
flows involving energy transfer [Chen and Adeniji-Fashola, 1987]. 


4. PARTICLE-TURBULENCE INTERACTION 


A very important aspect of the modeling of turbulent fluid-particle flows 
is the particle-turbulence interaction problem. Turbulence kinetic energy 
extracted from the mean flow kinetic energy of the continuous phase is partly 
dissipated by the smallest eddies and partly imparted to the particles thus 
enhancing the dispersion of the particulate phase. This ’two-way coupling’ 
referred to earlier - modulation of the kinetic energy of turbulence by the 
particles and enhanced dispersion of the particles by the turbulence will now 
be discussed in a little more detail. It is pertinent to point out at this 
point that the turbulent dispersion phenomenon is primarily responsible for the 
considerable enhancement in mixing observed for turbulent flows when compared 
with laminar flows. 


TURBULENT DISPERSION 

The turbulent dispersion phenomenon is very closely related to the inter- 
action between individual particles and turbulent eddies. A particle normally 
interacts with a series of eddies as it moves through the fluid. The particle 
trajectory scheme attempts to simulate this interaction by tracking each repre- 
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sentative computational particle through a succession of turbulent eddies con- 
tained within the domain of interest. Figure 1 is a schematic illustration of 

this interaction between particle and eddies and in relation to the computa- 
tional cells. As discussed by Gosman and Ioannides [1981], a particle interacts 
with a given eddy for a period of time which is the minimum between an esti- 
mated particle transit time within the eddy, t tr and an eddy lifetime, t e . The 

particle transit time is obtained as the solution of the linearized equation of 
motion of the particle while the Lagrangian time scale of the turbulent eddy is 

obtained from length and velocity scales of the turbulence which are extracted 

from a k- turbulence model. Thus, 


t int = Min[t e ,t (r ] 


(4.1) 


where 


ttr = -t* In [ 1 . 0-1 e /t*| Uj-vjl] 

and 


(4.2) 


4- , e/(2k/3 ) 1/2 


(4.3) 


The eddy length macroscale, l e is defined in terms of the kinetic energy 
of the turbulence, k and its dissipation rate, £ as 



r 3/4 3/2 
k Vs 


(4.4) 


In a stochastic formulation of the particle trajectory scheme which is the 
case in the present study, the fluctuating component of the fluid velocity, u’, 
is obtained from a Gaussian distribution of values having a zero mean and a 
standard deviation, 0- given by 


G ii - (2k/3)1/2 
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The Gaussian distribution is, however, not expected to be appropriate, in gen- 
eral, for describing non-homogeneous, non-turbulent flows. 


TURBULENCE MODULATION 

The presence of particles, even in very small concentrations, has the 
effect of modulating the turbulence intensity, the direction of modulation 
being influenced by the mean particle size and the level of modulation by the 
particle loading. This turbulence modulation effect was observed experimentally 
by Moderrass et al. [1984] and Tsuji et al. [1984] and attempts to mathemati- 
cally characterize the phenomenon include those of Al-Taweel and Landau [1977] 
and Chen and Wood [1985]. The interphase interaction force terms between par- 
ticles and the continuous phase are reflected as extra dissipation terms in the 
modeled equations for k and 6 when the former are included in the derivation of 
the field equations for the latter. The earlier attempts to implement these 
turbulence modulation models have been mostly within a two-fluid formulation in 
which the two phases are described as two interpenetrating continua viewed from 
an Eulerian framework. Equations (4.6) and (4.7) from Chen and Wood [1986] show 
the extra dissipation terms due to the turbulence modulation effect of the par- 
ticles for such a two- fluid formulation: 


a 


(«i *) = 



2k Pn 

— — 5- [1-exp (-0.5 t.£/k)] 
P 

K (TH2) 


(1+0.15R* p °- 687 H U l 

(TH1) 


-V,) 


(4.6) 


a 

dlj 


(U, e) = 


a 

ax, 


v, de e . K p e 

( i^- k lc ' p ^ £, -vr 

(TH3) 


(4.7) 


The term TH1 in equation (4.6) is the turbulence modulation term due to 
the mean slip while the terms TH2 and TH3 are due to the particle slip velocity 
at the fluctuating level. The model is valid for the situation 


t e > > l K » where t^ = (V/fc )*/ 2 


(4.8) 


is the Kolmogorov time scale. The model described above has been incorporated 
into the particle trajectory scheme of the present study. 


5. NUMERICAL SCHEME 


The set of governing differential equations describing the evolution of 
confined turbulent fluid particle flows cannot, in general, be solved analyti- 
cally thus requiring the adoption of a numerical procedure. For the continuous 
phase, the governing Eulerian equation set is solved using the SIMPLE algorithm 
of Patankar and Spalding [1972] and Patankar [1980], The overall scheme adopted 
for the solution of the governing equations is similar to that suggested by 
Crowe et al. [1977] and illustrated in Figure 2. An alternative scheme more suited to 
time-dependent flows was later presented by Dukowicz [1980] and further developed 
by Cloutman et al. [1982] and Amsden et al. [1985], 

First, the clean fluid flow field is obtained by solving 
the continuous phase governing equations. This is done using a staggered grid 
distribution in which velocity cells are centered about the edges of the scalar 
cells. Next, particle trajectories are computed for a predetermined number of 
representative particles such that a statistically stationary solution is 
obtained for the overall particle flow field. The particle trajectories, and 
temperature history where appropriate, are obtained by solving for the particle 
the non-linear ordinary differential equations of motion and the energy equa- 
tion subject to the currently existing continuous fluid flow and temperature 
fields. A fourth order Runge-Kutta algorithm is used for this purpose. During 
the calculation of a particle’s trajectory and temperature history, the sources 
of momentum, energy, kinetic energy of turbulence and its dissipation rate, all 
due to the particle motion, are accumulated for each computational cell trav- 
ersed. The form of the source terms have already been presented elsewhere [Ade- 
mji-Fashola and Chen, 1987] and so will not be repeated here. These source 
terms are then used in the next global iteration on the continuous phase field 
equations until convergence is attained. It was found that source term relax- 
ation was required to achieve stability of the global iteration scheme for some 
of the example flow problems studied. 


PARTICLE SOURCE FIELD CONTINUITY 

A necessary condition to obtain a globally converged solution is to ensure 

the continuity of the source fields as was also pointed out by Durst et al. 

[1984], In order to ensure compliance with this important requirement, it is 

necessary to ensure the computation of source terms for each cell traversed by 
each computational particle through a judicious choice of the particle integra- 
tion time step as well as have particles start from as many locations as is 
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practicable within the relevant portion of the inlet plane. In the present 
study particles are uniformly distributed in physical space at the appro- 
priate’ portion of the inlet plane of the computational domain in contrast to 
the scheme of Durst et al. [1984], in which particles are introduced only at 
grid nodes. The smooth profiles they obtained are very likely to be a 
consequence of the deterministic nature of the particle trajectories used in 

their study. 


INTEGRATION TIME STEP 


The choice of appropriate time steps for the integration of the particle 
equations of motion is very vital to obtaining a globally converged solution 
and smooth averaged particle flow fields. For the complex confined turbulent 
fluid-particle flow problems in general, some of the relevant time scales 
include the Lagrangian or macro time scale (eddy lifetime) of the turbulence, 
t • the Kolmogorov or the micro (dissipation) time scale of the turbulence, t K , 
the particle relaxation time, t*; the particle residence time within a computa- 
tional cell or the whole computational domain t R . Also relevant to the stochastic 
determination of the particle turbulent intensity are the particle transit time within an 
eddy t tP and the particle eddy interaction time, t int . The integration time step is 
selected to ensure adequate resolution with regard to the trajectory and temperature 
evolution while ensuring computational efficiency by avoiding unnecessarily small time 

steps. 


In the present study, a variable integration time step scheme was devised. 
An upper bound on the time step through any computational cell was imposed 

based on an estimated particle residence time for that cell and with the par- 
ticle being constrained to undergo about four integration steps within t e 
cell. Without this restriction, the possibility of a particle overshooting one 
or more cells, possibly due to a sudden reduction in cell dimensions in a 

non-uniform grid domain, exists. Such a situation will result in a failure to 
compute the relevant source term contributions for a cell that was actually 

traversed by the particle. The consequence will be a lack of smoothness m t e 
particle source distribution and, possibly, divergence of the global iter- 

ations. 


Also, for the reason of ensuring a smooth evolution of the particle tra- 
jectory and temperature history, a further restriction on the integration time 
step. At <t* , is made. The particle-eddy interaction time is determined and 
controlled independently of the integration time step. 
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PARTICLE AVERAGED PARAMETERS FROM PARTICLE TRAJECTORY STATISTICS 

One of the problems associated with the use of the Lagrangian particle 
trajectory approach, highlighted by Smoot and Smith [1985], is the difficulty 
of extracting smooth mean particle flow and temperature fields from the statis- 
tics of trajectories and temperature histories obtained for representative com- 
putational particles. In the present study, the fluid properties utilized in 
the particle trajectory and temperature history calculations are the linearly 
interpolated values in which the four nearest neighbors regarding the par- 
ticle s current location are used, resulting in second order accuracy [Sirig- 
mano, 1983], The details of the extraction of particle mean flow and tempera- 
ture fields information from the particle trajectory and temperature history 
statistics are available in Adeniji-Fashola et al. [1988]. 


BOUNDARY CONDITIONS 


The definition of a fluid flow problem becomes unique through the speci- 
fication of the boundary conditions after the governing differential equations 
are outlined and an appropriate closure of these equations is effected. The 
example flow problems investigated in the present study include vertical pipe 
upflow and horizontal recirculation chamber flow. However, rather than define 
the boundary conditions specific to each flow problem separately, the more 
efficient approach of defining generic boundary condition types is adopted. It 
then becomes a straightforward exercise to construct the boundary conditions 
for these and other specific flow situations of interest. 

Inlet Plane: 

The specification of the inlet plane boundary conditions for fluid flow 
problems is very important, as was discussed by Sturgess et al. [1983] and 
Westphal and Johnston [1984], since this influences significantly the 

subsequent evolution of the flow, especially in the case of parabolic flows for 
which the inlet plane conditions constitute the initial conditions for the 

solution of the governing differential equations. 

In order to correctly simulate a given fluid flow experiment numerically, 

the ideal specifications for the inlet flow variables are the actually measured 
values. The complete set of measured inlet flow variables is, however, hardly 

ever available. In the absence of such detailed experimental information, uni- 
form profiles are commonly specified for the axial velocity and temperature 
profiles of the continuous phase flow at the inlet plane. The turbulent kinetic 

energy is usually assumed to be a percentage, between 3 and 20%, of the inlet 

flow mean kinetic energy. The kinetic energy dissipation rate at the inlet is 
then obtained as 


64 



(5.1) 


E = ( CjV4 k 3/2 ) / 1(j 

where 1^, the dissipation length scale, is specified as a fraction of the 

characteristic length scale at the inlet. 

For the particle trajectory and temperature history calculations, the 
initial velocity and temperature slip values relevant to the particular flow 
problems are employed in setting the required inlet conditions. 

Exit Plane: 

At the exit plane, the usual boundary condition imposed for any flow var- 
iable, <{> , is 3^/ 9n = 0, where n is the normal to the exit plane. This con- 
dition is generally valid if the extent of the computational domain in the 
primary flow direction is sufficient to ensure fully-developed flow conditions 

for internal flows or self-similarity for jet flows at the exit plane. Particle 
trajectory and temperature history computations are discontinued for a compu- 

tational particle once the particle exits from the computational domain through 
the exit plane or any other open boundary. 

Solid Boundary: 

The conventional wall functions approach is used to impose wall boundary 

conditions on the velocity and temperature as well as the turbulence kinetic 

energy and its dissipation rate. The presence of particles in a fluid flow has 

been experimentally observed to influence the boundary layer [Kramer and Depew, 
1972] and, as a consequence, the nature of the wall function which is normally 

used to connect the actual value of a given variable at the wall to the value 
at the wall-adjacent grid node. During their trajectories, particles that reach 
the wall either adhere to it as observed in particle erosion problems [Dosanjh 
and Humphrey, 1984], or collide with the wall and get Reflected” back into the 
flow domain, usually with an accompanying loss of energy and momentum to the 
wall. In addition, the high level of shear in the wall vicinity coupled with a 
particle velocity slip introduces an additional transverse force on the par- 
ticle which further modifies its subsequent trajectory and behavior in the 

near-wall region. These effects have not been included in the present study, in 
which perfectly reflecting boundary conditions have been adopted for the parti- 
cle-wall interaction, but will be the subject of a future study. 

Other generic boundary condition types include the symmetry axis, for 
which d<$/dn = 0, where in this case, n is the normal to the symmetry axis, 
and the open boundary condition which has been used by Leschziner and Rodi 
[1984], Dosanjh and Humphrey [1984], Amano and Brandt [1984] and Chen and Ade- 
niji-Fashola [1987] for modeling parabolic flows of free jets and wall jets using 
elliptic formulations. These are described in greater detail by Adeniji-Fashola et al. 
[1988]. 
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6. EXAMPLE FLOWS 


In order to illustrate the multiple-realization particle trajectory mode- 
ling scheme for confined turbulent fluid-particle flows described above, two 
example flow problems - vertical pipe upflow and horizontal coaxial jet flow in 
a recirculation chamber with a particle-laden central jet and a clean annular 
jet are examined. 1500 computational particles were found to be adequate in 
each example for obtaining statistically stationary solutions. Typically, 
global under-relaxation values of 0.50 were found adequate to ensure the sta- 
bility of global iterations of which between five and seven were required to 
obtain globally converged solutions. The results obtained for the numerical 
simulation of these flows will now be discussed. 


VERTICAL PIPE UPFLOW 

The experimental data which served as the basis for this example numeri- 
cal simulation are those of Tsuji et al. [1984] for the upflow of a particle- 
laden stream in a straight vertical pipe. The experimental flow within the test 
section is considered to be fully-developed after going through a riser that is 
167.5 diameters long. 

A 50 X 23 uniform grid distribution was used to discretize the computational 
domain which had an axial extent of 60 pipe diameters. Figure 3 shows both the 
experimental data and the numerical predictions of the radial profile of the slip in 
the axial velocity between the air and the particulate phase. The mean particle size 
and loading ratio are 200 pm and 1.0 respectively. The air velocity is slightly 
overpredicted in the 0.2R - 0.8R range where R is the pipe radius. However, the 
prediction accuracy is considered to be good for such a complex system. The radial 
profile of the axial velocity of the solid phase is particularly well predicted. The 

location of the cross-over in sign of the slip between the phases is predicted to be 
closer to the wall, less than 0.1R from the wall, than the 0.2R from the wall that 
was experimentally observed. 

A similar picture obtained for the higher loading ratio of 2.1 is pre- 
sented in Figure 4. The level of accuracy of the predictions is similar to that 

of the 1.0 loading ratio case. However, it is the air velocity profile that is 

better predicted in this case. The solid phase axial velocity is considerably 

underpredicted in the inner 60 percent of the wall region. 

As pointed out earlier, the particulate phase has the effect of modulating 
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the level of the turbulence intensity. For smaller particle sizes this results 
in a decrease in the kinetic energy of turbulence. The experimentally observed 
and numerically predicted turbulence modulation effect for a loading ratio of 
3.2 are illustrated in Figure 5. The solid line in the figure shows the pre- 
dicted radial profile of the turbulence intensity for the corresponding clean 
flow. The predicted level of turbulence intensity is considerably higher than 
the level observed from experiment. Also, while a greater modulation effect was 
observed closer to the wall region, the predictions show a reversal in which 
the greater level of modulation is located closer to the pipe centerline. The 
imposed wall boundary conditions and wall functions in the numerical scheme are 
probably responsible for the suppression of the modulation effect in the 

near-wall region. 

The development in the axial direction of the streamwise velocity of the 
particulate phase for an inlet velocity slip ratio of 0.10 is shown in the 

contour plot of Figure 6a and a corresponding surface plot in Figure 6b The 

ability of the particle trajectory scheme to effectively handle extreme leves 
of velocity slip was tested by imposing an axial slip velocity of 0.10 at t e 

pipe inlet plane. The figures indicate that a fully developed state was 
attained in the 60D extent of the computational domain. 


HORIZONTAL COAXIAL JET FLOW IN RECIRCULATION CHAMBER 

The experimental data of Leavitt [1980] serve as the basis for the numer- 
ical simulation of this example. The actual geometry studied is illustrated in 
the schematic of Figure 7. The primary jet air velocity at inlet is 33 m/s 

while the corresponding secondary jet air velocity is 42 m/s. Coal particles of 
a mass mean diameter of 43l*m were used to uniformly seed the primary jet and 
the particle loading ratio is 1.50. The estimated turbulence intensity levels 
at the inlet are 15 and 18% for the primary and secondary jets respectively. 
The primary and secondary jet diameters at inlet are 0.0255m and 0.127m respec- 
tively while the chamber diameter is 0.206m. The axial extent of the recircula- 
tion chamber is 0.926m (36.3 primary jet diameters or 4.5 chamber diameters). 

A 41 X 41 non-uniform staggered grid distribution, shown in Figure 8, is 

used for the numerical study and the computational domain extended to 20D where 
D is the chamber diameter. The numerical prediction of the evolution of the 
axial velocity is shown in Figure 9. The corner recirculation zone is seen to 
extend to about 1.79D. No particles are predicted as reaching this recircula- 

tion zone and this is believed to be due to the high chamber-to-primary jet 
diameter ratio of 8.08 and the positive slope of the shear in the mixing layer 

between the primary and the secondary jets which will result in a slip-shear 
transverse force directed towards the centerline. Another interesting observa- 
tion is that the particle axial velocity starts lo lead that of the fluid from 
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about the 0.80D axial location and this continues to about the 7 43D axial 
location downstream of which all axial velocity slip disappears. Particles are 
seen to have dispersed to the outer extremities of the recirculation chamber by 

the time the 12.0D axial location is reached. However, it should be remembered 

that this is only a hypothetical situation since the actual experimental inves- 
tigation was limited to an axial extent of only 4.5D. 

Figure 10 shows the axial evolution of the turbulence intensity It is 
observed that up to about the 3.0D axial location, the turbulence intensity in 

the presence of particles (shown dotted) falls below that of the clean flow in 

the primary jet portion of the flow but is actually higher for the rest of the 

chamber in the radial direction. However, beyond the 3.0D axial location, the 
clean flow turbulence intensity uniformly lags the two-phase intensity at all 
radial locations for any given axial location. The kinetic energy of turbulence 
is essentially fully developed at the 5.15D axial location and only a radially 
uniform decrease in magnitude is observed for the rest of the flow in the axial 

direction. This is in contrast to the radial profile of the axial velocity 

which does not become fully developed for both phases until the 12.0 to 15 0 

diameter range is reached. 

The contour and surface plots of the particle axial velocity are shown in 
Figures 11a and lib. These have been normalized with respect to the secondary 

jet gas velocity at inlet. Since, in contrast to the two-fluid scheme, non-zero 

values of the particle velocity are not returned for computational cells not 
visited by any particle during the trajectory calculations, the zero-velocity 

surface in the plots of Figure lib also indicate the particle-deficient 

regions. 


The comparison of the limited experimental data available from Leavitt 
[1980] is currently being undertaken. 


7 - c oncluding remarks an d recommendations for ftjturf. work 


The numerical modeling of confined turbulent fluid-particle flows using 
the multiple-realization particle trajectory scheme has been presented. The 
performance of the numerical modeling scheme has been tested using data for the 
upward flow of a fluid-particle stream in a straight vertical pipe and for the 
horizontal coaxial jet flow in a large recirculation chamber for which the cen- 
tral jet is particle-laden. 


The multiple-realization particle trajectory turbulent flow modeling 
scheme ... 6 
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is believed to be a more elegant and accurate approach to the 
extension of single-particle hydrodynamics to dilute multi- particle 
systems; 

is better able to incorporate additional force terms such as lift, 
virtual mass and Bassett history terms in the particle equation of 
motion as appropriate; 


needs further investigation in order to improve its computational 
efficiency and so reduce its huge CPU time requirements; 

needs to have the particle-turbulence and particle-wall interactions 
further investigated to improve prediction accuracy. 
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PARTICLE DISPERSION MODELS AND DRAG COEFFICIENTS 
FOR PARTICLES IN TURBULENT FLOWS* 

by 

C.T. Crowe, J.N. Chung and T.R. Troutt 
Department of Mechanical and Materials Engineering 
Washington State University 
Pullman, WA 99164-2920 

INTRODUCTION 

The dispersion of particles in turbulence is fundamental to a variety of mass apd en- 
ergy transfer processes. The dispersion of particles in jets is important to the combustion 
process and the design of propulsion systems. The separation of particles by electrostatic 
precipitation is important in many applications from the power industry to clean-room 
technology. Understanding the basic phenomena underlying particle transport in turbu- 
lence and establishing viable models is important to the development of new technologies 
for advanced propulsion systems. 

This paper reviews some the concepts underlying particle dispersion due to turbu- 
lence. The paper addresses the traditional approaches to particle dispersion in homoge- 
neous, stationary turbulent fields and reviews recent work on particle dispersion in large 
scale turbulent structures. The paper also reviews the state of knowledge on particle drag 
coefficients in turbulent gas-particle flows. 

MECHANISMS FOR PARTICLE DISPERSION 
Basically two mechanisms have been used to physically model particle dispersion in 
turbulence. In homogeneous turbulent flows, the most common approach is to regard 
dispersion as a stochastic process. On the other hand, dispersion in large scale turbulent 
structures appears to be more influenced by the ordered motion. Both approaches are 
treated separately. 

Particle Dispersion in Homogeneous, Stationary Turbulence 

The traditional approach to treating particle dispersion is turbulence is to regard 
the process as a gradient diffusion (Fickian) process in which the diffusional velocities 
are proportional to the concentration gradient, the constant of proportionality being the 
diffusion coefficient. 

dc 1 \ 

= ~ D tor '> 

The earliest work that related the diffusion coefficient to properties in a homogeneous, 
stationary and isotropic flow was that of Taylor (1921) who provided the following rela- 
tionship, 

D = d*T L 2) 
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where u 1 ^ is the mean square of the fluctuating velocity in the direction transverse to 
the main flow direction and Tj, is the lagrangian integral time scale. Although Taylor’s 
analysis was done for a fluid point, Synder and Lumley (1971) showed it was applicable 
to a particle provided the lagrangian time scale corresponded to the particle trajectory. 

The earliest study of particle motion in a turbulent field was reported in a PhD thesis 
by Tchen (1947) who integrated the Basset, Boussinesq, Oseen equation for a particle in 
a homogeneous, stationary turbulent field. He assumed that the particle remained always 
in the same turbulent eddy. By so doing the long-time diffusion of the particle was equal 
to that of a fluid particle. 

Many researchers after Tchen have strived to improve Tchen’s model by relaxing the 
assumptions made by Tchen. Peskin (Soo, 1967) solved a nonlinear stochastic equation 
for the motion of a particle which did not deviate far from the initial coincident turbulent 
eddy. Unlike Tchen’s analysis, Peskin assumed that only Stokes drag acted on the particle. 
He predicted that the diffusivity decreased uniformly with an inertial parameter which 
related the aerodynamic response time of the particle to the time scale for turbulence. 
The physical argument underlying this result was the larger the aerodynamic response 
time of the particle is compared to the eddy life time, the less a particle would respond to 
the unsteady turbulent field. Hinze (1972) used similar time-scale arguments for particle 
dispersion and concluded that if the particle density ratio is large, only those particles less 
than one tenth of the dissipation length scale will respond to the turbulent fluctuations. 

A more general analysis of particle dispersion in homogeneous, isotropic stationary 
turbulence has been reported by Reeks (1977). He assumed a linear drag law and body 
force acted on the particle and obtained an expression for diffusion which depended on 
time, particle aerodynamic response time and the correlation function for the velocity 
field. He then utilized Phythian’s model (1975) for the turbulent energy spectrum and 
predicted a particle diffusion coefficient. His results show that the diffusion coefficient for 
particles with no body force increases with increasing time and, at long times, approaches 
an asymptotic value. This long-time diffusion coefficient increases with increasing aero- 
dynamic response time and can exceed that of a fluid particle. This result differed from 
Peskin’s model. The reason underlying the trend predicted by Reeks is the fact that 
the diffusion is jointly dependent on the mean square of the particle velocity fluctuations 
and lagrangian integral time scale as shown by equation 2. Although the amplitude of 
the fluctuation velocity of the heavier particle is reduced, the lagrangian time scale is 
increased proportionately more giving rise to an increased diffusion coefficient. 

Another factor controlling particle dispersion in homogeneous, stationary turbulence 
is the “crossing trajectory” effect first identified by Yudine (1959). If the mean velocity of 
the particles is different from that of the fluid, such as particles dropping at their terminal 
velocity through a turbulent field, the particles remain less time in a given eddy. The 
reduction in fluid-particle interaction time reduces the particle diffusion coefficient. Reeks 
also predicted a decrease in long-time diffusion coefficient by including a body force in 
the particle motion equation. 

Thus there are two primary factors controlling particle dispersion in homogeneous 
turbulence, the inertial effect and the crossing trajectory effect. Experimentally it has 
been difficult to separate these effects since a heavy (not a fluid) particle will have both 
inertial and crossing trajectory effects. The most convincing experimental evidence that 
separates the two effects have been provided by Well and Stock (1983). They suspended 
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glass beads by Coulomb forces in a horizontal flow of near-homogeneous, grid-generated 
turbulence. The crossing trajectory effect was controlled by adjusting the field strength. 
Their data show that the inertial effects on particle diffusion coefficient are small compared 
the crossing trajectory effect. The small increase in diffusion with increasing inertial effects 
predicted by Reeks was not discernible because of the scatter in the experimental data. 

Future developments in the analysis of stochastic turbulent flows will address depar- 
tures from homogeneity and isotropy. Reeks (1981) has initiated work in this direction. 

Particle Dispersion in Large Scale } Turbulent Structures 

Large scale turbulent structures are encountered in flows generated by a large ve- 
locity gradient such as free shear layers and jets. Under these circumstances, large scale 
turbulent structures are formed which grow and pair with time. These structures were 
first identified by Brown and Roshko (1974) in flow visualization studies of mixing layers. 
A typical photograph of a large scale structure is shown in figure 1 . These turbulent 
flow fields are inhomogeneous, non-stationary and anisotropic but represent important, 
practical problems in industrial applications such as combustion systems. 

Particle dispersion in turbulent flow characterized by large scale structures is mech- 
anistically different than that in homogeneous flows. The particle motion is controlled by 
the moving structures and not by the fine scale turbulence. Thus, the dispersion process 
cannot be regarded as gradient transport. 

A conceptual model for particle dispersion ij> large scale structures is the entrapment 
of particles in the structure and the subsequent centrifuging of the particles beyond the 
structures. This concept was first suggested by Singamsetti (1966) who observed, ex- 
perimentally, that particles in a submerged jet could disperse more quickly than a fluid 
particle. The same trends have been observed by Lilly (1973), Householder 1968), Laats 
and Frishman (1970) and Subramanian and Ganesh (1984) in the experimental study of 
particle and droplet laden free jets. Lilly attributed his results to an increase lagrangian 
time scale but Yuu et al. (1978) claim Lilly’s results were a manifestation of his experi- 
mental set-up. Laats and Frishman noticed this trend only in the early portion of the jet 
formation and surmised that it was due to a Magnus effect. Goldschmidt et al. (1972), in 
reviewing Householder’s data, mentioned the possibility of particle centrifuging by large 
scale structures but concluded that the mechanism was not viable because it did not 
explain the observed trends in centrifuging bubbles. 

Yule (1981) gave some credence to the mechanism when he observed droplets in jet 

flows being centrifuged toward the outer flow. 

Crowe et al. (1985) report an effort to quantify those conditions under which the large 
scale structures are responsible for dispersing heavy particles beyond fluid particles. They 
proposed a time scaling argument similar to that used by Hinze (1972) for homogeneous 
turbulence. The aerodynamic response time of a particle is the time required for a particle 
released from rest in a uniform flow to accelerate to 63% of the flow velocity. For Stokes 

flow it is 


pd 2 

TA = 18 p 


3 ) 


where p is the fluid density, d is the particle diameter and p is the dynamic viscosity of 
the fluid. It is a measure of the responsiveness of a particle to changes in flow velocity. 
The characteristic time of the flow is given by 


f 
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Figure 1. Shadowgraph Visualization of Large Scale Turbulent Structures in Plane 
Mixing Layer (Brown and Roshko, 1974) 
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where S is the mixing layer thickness and A U is the velocity difference across the layer. 
Thus the scaling factor is 


T A pd? AU 

Tp 18 p.8 
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A schematic diagram showing the effect of Stokes number on particle dispersion 
is shown in figure 2. For small Stokes numbers, the particles are in near equilibrium 
with the conveying fluid and the particles will disperse as a fluid particle. At large Stokes 
numbers, the particles have insufficient time to respond to the structures and they disperse 
less than the fluid particle. However at intermediate Reynolds numbers, the particles are 
entrapped by the rotating structure and are centrifuged beyond the structures giving rise 
to a dispersion exceeding that of a fluid particle. 

A preliminary numerical study (Crowe et al. 1985) using a simple vortex sheet model 
proposed by Stuart (1967) showed that particles could be centrifuged beyond the vortex 
at Stokes number between 0.1 and 1 as shown in figure 3. A subsequent study by Gore 
et al (1985) using pseudospectral direct simulation for modeling the vortex structures 
showed the same trend lending support to the model. 

Recent experimental studies have shown the importance of large scale structures in 
the turbulent dispersion process. Kamalu et al. (1987) at Washington State University 
have reported on particle dispersion studies in horizontally oriented plane mixing layer. 
Both naturally evolving and subharmonically forced mixing layers were studied. The 
forced mixing layer was generated by a sound source, the effect of which was more ordered 
vortex structures. Particles, 40 microns in diameter, were released from a source upstream 
of the layer and photographed. A photograph of the particles in a subharmonically forced 
mixing layer is shown in figure 4. Shown in the same figure are streaklines generated by 
smoke with no particles in the flow. One notes the absence of particles in the vortex cores 
and the accumulation of particles near the edge of the vortex structures. 

Laser Doppler measurements of particle velocities by Wen et al. (1987) lend more 
quantitative support on the role of large scale structures in particle dispersion. The 
measured lateral particle and fluid velocities taken in the same facility used by Kamalu 
et al. are shown in figure 5. The fluid velocities show the expected trend, negative on 
the top (high speed side) and positive on the bottom (low speed side) which indicates 
a motion of the fluid towards the mixing region. The particle velocities, on the other 
hand, are positive on the top and negative on the bottom indicating motion away from 

the mixing layer. . 

The importance of large scale structures in particle dispersion has also been verified 

in recent experimental studies by Kobayashi et al (1987) and Lazaro and Lasheras (1987). 

The importance of large scale structures in turbulent dispersion of particles has been 
established. It represents a demarcation from particle dispersion in homogeneous turbu- 
lence because it is more deterministic than stochastic. Therefore it is not reasonable to 
model particle dispersion as a gradient transport process. 
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Figure 2. Schematic Diagram Illustrating the Effect of Large Scale Structures on 
Particle Dispersion 
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Figure 3. Predicted Particle/Fluid Dispersion Ratio as Function of Stokes Number 
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Figure 4. Photographs of Fluid Streaklines and Particles in Large Scale Turbulent 
Structures (Kamalu et a/., 1987) 







NUMERICAL MODELS FOR PARTICLE DISPERSION 

Numerical models for fluid-particle flow can be divided into two categories; two-fluid 
models and trajectory models. A review of these approaches has been provided by Crowe 
(1982). In the two fluid model, the particulate phase is regarded as another fluid. In 
the trajectory model, the particle field is established by integrating particle trajectories 
through the field. 

Gradient Transport Models 

It is natural in the two fluid model to treat particle dispersion as a gradient transport 
process and assign a diffusion coefficient to represent particle dispersion due to turbulence. 
Elghobashi et al. (1983) have considered in detail the two-fluid model with the two- 
equation model for turbulence and applied it to free jets. They recognize that the gradient 
transport assumption is valid only if the energy containing eddies are much smaller than 
the length scale for the transport gradient. They suggest a correction term to Ficks law 
that represents a convective flux due to flow inhomogeneity and which disappears for 
homogeneous turbulence. The diffusion coefficient is related to the kinematic viscosity 
through an effective Schmidt number, the value for which is not provided in the paper. 
The model requires three additional empirical constants above those required for the k — e 
model in single phase flows. They apply their model to the prediction flow properties of a 
jet studied experimentally by Modarress et al. (1984) and claim good agreement between 
measurements and predictions. 

Chen and Wood (1985) also use the two-fluid approach to model a jet. They as- 
sume that the particle and fluid phase have the same average velocity and justify this 
assumption on time scale arguments. They also use a gradient transport assumption for 
particle dispersion due to turbulence and chose an effective Schmidt number of 0.7 which 
corresponds to the value for diffusion of a passive scalar in a round turbulent jet. Sub- 
ramanian and Ganesh (1984) measure an effective Schmidt number of 0.47 for the same 
configuration. Chen and Wood applied their model to experimental studies reported 
by Modarressef al. (1984) and Wood et al. (1984) and noted good agreement between 
predictions and measurements. 

The real difficulty in using the gradient transport model for particle dispersion in 
connection with the two-fluid model is encountered in wall-dominated flows. Here one 
has to establish boundary conditions for the particle concentration and velocities. The 
concentration profile will depend on the interaction of particles with the wall. For example, 
if the particles stick to the wall a different boundary condition must be used than if the 
particles rebound elastically from the wall. For inelastic collisions, another assumption 
must be used. 

In addition the particle velocity component parallel to the wall is not zero as in a 
single phase continuum flows. Chen (1986) utilizes concepts from rarefied gasdynamics 
and calculates a slip velocity which depends on a fluid-particle interaction length. The 
normal component of velocity is set equal to zero although this would not be true for 
an inelastic collision. Other ramifications of the two-fluid versus trajectory models are 
discussed by Crowe (1986). 

The advantage of the gradient diffusion model for particle dispersion is that it can 
be accommodated directly into the two fluid model. It also does not require excessive 
computational times. The difficulty is the selection of appropriate Schmidt numbers 
and other empirical parameters for a given application. The relative advantages and 
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disadvantages of the eulerian and trajectory approaches are discussed by Durst et al 
(1984). 

Monte Carlo Methods 

The method most natural to predicting particle dispersion using the trajectory ap- 
proach is the Monte Carlo method. By this method the turbulent field is represented 
with a random number generator. The basic idea was first proposed by Hutchinson et al. 
(1971) and was subsequently used by Gosman and Ioannides (1981) in conjunction with 
the k — e turbulence model for sprays. The turbulent fluctuational velocity is selected 
from a random number generator with a variance proportional to the turbulence energy. 
The particle motion equation is integrated with this velocity field until it passes from the 
eddy. The particle-eddy interaction time is established by the characteristic life time of 
the eddy or by the time for the particle to pass through the eddy. The dissipation length 
and time scales are chosen as the characteristic size and time and are given by 

L e = Cj /4 A 3 ^/« 6a) 

T e = L e /(2k/2) 1/2 66) 

where k is the turbulence energy, e is the dissipation rate and C p is an empirical constant 
arising from the k — e model. The time to pass through an eddy is approximated by 



where U p is the particle velocity and U g is the mean gas velocity. The interaction time 
is the minimum of the eddy life time and the passage time. If the passage time is small 
compared to the eddy life time then the crossing trajectory effect is important. 

Gosman and Ioannides applied this model as a test case to the experimental study 
done by Snyder and Lumley (1971) who measured the dispersion of a series of particle 
types in grid-generated turbulence produced by a vertically oriented wind tunnel. Synder 
and Lumley found that the heavier copper particles dispersed less than the lighter hollow 
glass beads. From the current state of knowledge, it is accepted that this trend is due to 
the crossing trajectory effect. Gosman and IonEfddies report good agreement between the 
their predictions and Synder and Lumley’s results even though their droplet equations 
do not contain a body force term due to gravity. This was probably an omission in the 
paper. 

The Monte-Carlo technique was used by Chen and Crowe (1984) to model particle 
dispersion measurements in fully developed pipe flow reported by Arnason and Stock 
(1984). As with Gosman and Ioannides, they found that the technique worked well for 
the near isotropic, homogeneous field in Synder and Lumleys experiments. However, it 
was necessary to tune the model by changing C p to achieve the best fit. Applying the 
model to the pipe flow experiments yielded very poor agreement with the experimental 
results as shown in figure 6. The model predicted that the larger particles would disperse 
less than the small particles due to the crossing trajectory effect but the opposite trend 
was found experimentally. Chen and Crowe rationalized that the turbulence model was 
too crude for the complex turbulent flow in a pipe. One shortcoming of the model is the 
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lack of information on the lagrangian length scale which should be used for L e defined 
above. The above equations provide, at best, an estimate of the lagrangian length scale. 

Faeth et al. (Solomon et al . , 1983; Solomon et al, 1985) have used the Monte-Carlo 
model extensively with their model for particle and droplet laden jet flows and show good 
agreement with experimental results. The model is calibrated to fit the analytic results of 
Hinze (1975) for diffusion of fluid particles from a point source in homogeneous, isotropic 
turbulence. 

The general utility of the Monte-Carlo method for particle and droplet dispersion 
remains to be established. The Monte Carlo method is attractive because of the mini- 
mal empiricism needed to model the flow (provided the stochastic representation of the 
turbulence is reliable) and the simplicity of handling boundary conditions. The primary 
problem is the large number of trajectories needed to establish a stationary average in a 
computational cell. 

Another dispersion model which attempts to capture the desirable qualities of the 
gradient transport and trajectory approaches is the “hybrid” model first proposed by 
Jurewicz (1976) and subsequently used by others. In this approach the trajectory model 
is first used to calculate particle concentrations in each cell. Then, a diffusional velocity 
is added to the particle velocity which is proportional to the concentration gradient and 
diffusion coefficient. Of course, this model requires selection of a diffusion coefficient. 

Nonstationary, nonhomogeneous models 

The gradient transport models are inadequate to model particle dispersion by large 
scale structures. Since particle dispersion appears to be controlled by the large scale 
motion, the numerical model must represent the essence of these flows. Work in this area 
is just beginning to appear in the literature. 

Chein and Chung (1987) report a report a numerical study of particle dispersion in 
vortex pairs modeled using discrete vorticies. They found that particles with intermediate 
Stokes numbers (0.5 to 5) are dispersed more than the fluid particles while at larger Stokes 
numbers the heavy particles disperse less than the fluid particle. 

Chien and Chung (1988) also used the discrete vortex method for generating a time 
dependent two-dimensional mixing layer. Particles, released in this flow field show the same 
general trend as shown by figure 7. At low Stokes numbers, the particle? follow the fluid 
and disperse as a fluid particle. As the Stokes number is increased the particles disperse 
more and near a Stokes number of unity the vortex core is almost void of particles as they 
are centrifuged out. In this regime the particle dispers^more than the fluid particle. With 
further increase in Stokes mumber, the particles become unresponsive to the vortices and 
more in near rectilinear trajectories. The same trend is noted in a numerical study of jets 
by Chung and Troutt (1988). 

Future work in particle dispersion in large scale structures will witness more ad- 
vanced fluid mechanic models such as vortex models and pseudo-spectral methods as the 
computational capability is enhanced by future generation computers. 

PARTICLE DRAG COEFFICIENTS 

Fundamental to the development of numerical models for gas solids or gas-droplet 
flows is the particle or droplet equation of motion. In general, there are many forces 
acting on the disperse phase particle such as the virtual mass force, Basset force, pressure 
gradient force and the steady state drag force. The most widely accepted formulation for 
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Figure 6. Measured and Predicted Particle Dispersion in Fully Developed Duct 
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the equation of motion for low Reynolds number flow is that of Maxey and Riley (1983 
who derived the equation from basic principles. For gas-solids flows in which the material 
density of the particle is three orders of magnitude larger than the conveying phase, .the 
primary force is the steady state drag force which is quantified by the value of the drag 
coefficient, C D , and related to the steady state drag by. 

Fj) = — pApC d{U g — ~~ Up\ 

where p is the gas density, A p is the projected area and (U g - U P ) is the relative velocity 
vector between the fluid and the particle. 

There is a plethora of literature available on particle drag coefficient. Most of the 
data have been obtained for single particles or spheres mounted in an airstream. However 
in numerical model development, one is more interested in the drag coefficients of particle 
in a cloud. The particle drag data show significant discrepancies as shown in figure 

Ingebo (1956) published a NACA report on particle drag coefficient which he mea- 
sured by releasing solid particles in an airstream downstream of a grid. The particles 
were tracked by a rotating mirror camera and the velocity-distance data were reduced to 
obtain the acceleration and drag force. Ingebo found the drag coefficient was less than 
the standard value for a sphere and attributed the discrepancy to the acceleration of the 
particles. Crowe (1962) suggested that the low value could have been due to a critical 
Reynolds number effect created by the grid upstream of the particle injection ocatio • 
Arrowsmith (1973) suggested that the local air velocity in the cloud was less the 

tunnel speed affecting the calculation of the relative velocity. The discrepancy has yet 

be resolved. 

Hanson (1952) measured the deceleration of hexane droplets issuing from an atomizer 
into an air flow. The spray was photographed and the droplet deceleration reduced from 
the photographs. He assumed that the local gas velocity was constant throughout the 
chamber. Hanson’s results for C D are very low. Hanson attributed the low drag coefficient 
to evaporation but this explanation seems unlikely. 

Rudinger (1969) injected particles into a vertically oriented shock tube and passed a 
shock wave through the particle cloud. He used a rotating drum camera to record particle 

motion. The drag coefficients he reduced were significantly higher than the stand 

curve. Rudinger hypothesized that the turbulence generated by the particles created ig- 
zag motion which made the particles appear to have an higher “effective drag coefficient. 
However the same trend would have been noticed in Ingebo s results. 

Crowe (1962) reported measurements on the drag coefficient of burning gun powder. 
The burning powder was subjected to a shock wave in a vertical shock tube in a manner 
similar to Rundinger’s experiment. Particle motion was measured with a high speed 
camera. The data lie slightly above the standard curve but well below Rudinger s. 

Briffa (1981) has reported on the measurement of droplet drag in sprays. A water 
spray was photographed to yield a triple exposure of a droplet. The local air veloci y 
was measured by photographing the motion of Lycopodium dust The reduce g 

coefficients were smaller than the standard curve. Briffa attributed this tren o 
Basset term in the equation of motion but it is unlikely that the Basset term would be so 

predominant. 
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Tsuji et al. (1982) generated a stationary array of particles and measured the drag 
on one particle in the array using the pendulum method. Two configurations were tested; 
side by side particles and one particle in the wake of another. They found that the drag 
of the particle in the wake was reduced for separation distances of less than 10 diameters. 
The difference in drag coefficient does not explain the discrepancies observed in figure 8. 

Tsuji et al. (1984) reported on the LDV measurements of particle and air velocities 
of 200, 500 and 3000 micron particles in a vertical pipe. Small tracer particles were 
used to measure the gas flow velocity and the signals from the test particles and tracer 
particles were separated by a special signal discrimination device. The drag coefficients 
resulting from Tsuji et aV s experiments have been reduced by Lee (1987). The data 
fall below the standard drag curve but demonstrate significant scatter. Lee correlates the 
data with particle volume fraction, Froude number, Reynolds number (based on turbulent 
fluctuational velocity) and density ratio (particle to fluid material density ratio). By so 
doing, he was able to fit the data on a single curve. Still, extension of the empirical results 
to other conditions is tenuous because one would not anticipate that the aerodynamic drag 
would depend on the density ratio. 

Very recently, Fleckhaus et al. (1987) have reported measurements of particle ve- 
locities and concentrations in a jet with a two- dimensional LDA system. They also had 
tracer particles in the jet to obtain the gas-phase velocity. By fitting cubic splines to 
their velocity measurements, they were able to reduce particle accelerations. The drag 
coefficients were obtained by knowing the particle (glass beads) size and relative velocity. 
Their drag data lie above, but close to, the standard curve. 

There is a need to establish a valid drag coefficient for particles in a turbulent flow 
and to resolve the many discrepancies apparent in the data. It would be appropriate to 
repeat some of the earlier experiments using modern instrumentation to either validate 
the data or indicate the reason for the observed trends. Until more specific information 
is available, one is advised to use the standard drag curve for an isolated particle. 
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Recent work has demonstrated that the time dependent properties exhibited by concentrated 
suspensions of non-colloidal spheres when sheared in a conventional Couette viscometer may be 
explained in terms of shear-induced particle migrations (Gadala-Maria and Acrivos, 1980; 
Leighton and Acrivos, 1987b). These suspensions were observed to exhibit both a short-term 
increase in viscosity upon shearing immediately after loading into the Couette device and a 
subsequent long-term decrease after prolonged shearing, which were used to estimate the effective 
shear-induced diffusivity for concentrated suspensions both normal to the plane of shear and 
parallel to gradients in fluid velocity within the plane of shear (Leighton and Acrivos, 1987b). 

In this paper the experimental evidence for the existence of shear induced migration processes 
is reviewed and the mechanism proposed by Leighton and Acrivos (1987b) is described in detail. 
The proposed mechanism is shown to lead to the existence of an additional shear induced 
migration in the presence of gradients in shear stress such as would be found in Poiseuille flow, 
and which may be used to predict the amplitude of the observed short-term viscosity increase. 

The concentration and velocity profiles which result from such a migration are discussed in detail 
and are compared to the experimental observations of Kamis, Goldsmith and Mason (1966). 


1. Introduction 

Particle migrations across fluid streamlines in suspensions may result from a wide variety of 
mechanisms, ranging from Brownian type diffusive motions to the inertia induced drift 
mechanisms studied by Ho and Leal (1974) and others. In this paper we are concerned with 
shear-induced, particle migrations which have been observed to occur in concentrated suspensions 
of non-colloidal particles (particles sufficiently large that colloidal forces are unimportant at 
distances comparable to the particle diameter) and at sufficiently low Reynolds numbers that 
inertial forces may be neglected. Particle migrations under these conditions have been observed 
by a number of researchers. Early work by Kamis and Mason (1967) demonstrated that particles 
tend to accumulate behind an advancing meniscus in flow through a tube, and to be depleted 
behind a receding meniscus, the magnitude of the phenomenon being a strong function of the 
particle diameter/tube radius ratio and the concentration, suggesting that the particle migration was 
the result of some type of wall effect. Leighton (1985) demonstrated that if proper precautions 
were not taken, this effect can lead to serious errors in viscosity measurements in concentrated 
suspensions. The coefficient of shear-induced self-diffusion of spheres in a sheared suspension 
(i.e., the diffusion or dispersion arising from a random walk of particles in a sheared suspension 
analogous to Brownian diffusion) was examined by Eckstein, Bailey and Shapiro (1977) and later 
by Leighton and Acrivos (1987a). . 

The observations of particle migration of primary interest here were initially made by 
Gadala-Maria and Acrivos (1980) in suspensions of 40(im to 50|im diameter polystyrene spheres 
in silicone fluids. In the course of viscometric measurements of concentrated suspensions made 
with a conventional Couette viscometer, Gadala-Maria and Acrivos (1980) found that the 
suspension viscosity would decrease after prolonged shearing, and eventually reach a steady-state 
value which was as much as a factor of two below the initially observed value (cf. Figure 1). In 
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Fjcuhe 1 Relative viscosity of a <j> = 0.45 suspension as a function of the time that it had been 
sheared m the Couette device at y - 24 s' 1 . Polystyrene spheres, 40-50 pm in diameter in a mixture 
ot silicone oils (from Gadala-Maria 1979, figure 33). 


subsequent experiments, Leighton and Acrivos (1987b) demonstrated that the viscosity decrease 
was due to particle migration out of the sheared gap and into the reservoir by sealing the base of 
the Couette device with a layer of mercury and showing that the phenomenon disappeared. The 
rate of particle migration was found to be proportional to the shear rate and the square of the 
particle radius, and was successfully modelled by a one-dimensional diffusion process. The 
viscosity decrease was thus used to measure the effective diffusivity in the direction normal to the 
plane of shear for concentrated suspensions. The effective diffusivity was found to be a very 
strong function of concentration (cf. Figure 2) and to be much larger than the shear-induced 
coefficient of self-diffusion measured by Leighton and Acrivos (1987a). 

During the course of their experiments (also with polystyrene spheres in silicone oils) 
Leighton and Acrivos (1987b) observed that, upon first shearing the suspension in the Couette 
device, the viscosity would increase over a total strain of approximately 100, reaching a 
steady-state value before the subsequent long-term viscosity decrease. Since the long-term 
viscosity decrease only became significant after a strain of about 10 2 3 had elapsed, the two 
phenomena were well separated in time and could be investigated independently. As it is central 
to the present investigation into concentration distributions in Poiseuille flow, the short term 
viscosity decrease phenomenon is examined in more detail in the next section, which follows the 
development by Leighton and Acrivos (1987b). 


2. Short-Term Viscosity Increase 

The observed timescale appropriate to the initial viscosity increase phenomenon was found to 
be inversely proportional both to the shear rate y and to the square of the ratio of the particle 
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Figure 2 Diffusion coefficients calculated from the long-term viscosity decrease experiments: □, 
46 pm polystyrene in 0.639 mm gap; O. 46 pm polystyrene in 1.261 mm gap; A, 87 pm 
polystyrene in 1.261 mm gap. Shear rates were: filled symbols, 76 s -1 ; open symbols, 24 s -1 ; 
half-filled symbols, 7.6 s _1 . 


diameter to gap width ratio 3; thus this effect was explained in terms of a shear-induced migration 
of particles across the width of the Couette gap analogous to the long-term viscosity decrease 
phenomenon. This could only be the case, however, if the suspension flowing into the gap 
during the loading procedure acquired a concentration distribution across the gap that was 
different from the equilibrium profile corresponding to Couette flow. Then, upon shearing, the 
particle distribution would diffuse into that appropriate for the Couette flow, and thereby induce a 
change in the observed viscosity. The actual change in the viscosity results from the non-linear 
dependence of the observed viscosity on the concentration profile. For example, if the 
concentration profile is only slightly non-uniform and is assumed to be symmetric about the 
centerline of the gap, then the observed viscosity may be calculated to be: 


Ha* 

If 




i^iy 2 II 
If df ‘2 If 



(A<J>) 2 dy + 0(<( A<(>) 3 >) 


( 2 . 1 ) 


where the viscosity has been expanded in a Taylor series about <|) =m|>. In equation 2.1, |i is the 
viscosity that corresponds to a uniform concentration <fT, A<|> = is the deviation from the 
average concentration across the gap, y = 0 denotes the centerline of the gap, and y = ±b the 
walls. Note that the variation in viscosity is proportional to the average value of (A<J>) 2 across the 
gap for small fluctuations in concentration. The viscosity function in equation 2.1 may be 
calculated from the dependence of viscosity on concentration observed by Leighton and Acrivos 
(1987b): 
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an Eiler's equation where <|> m is the maximem particle concentration and [|i] is the intrinsic 
viscosity. The best fit values for the suspensions used in the experiments discussed here in the 
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range .3 < <J> < .5 were <}> m = 0.58 and [|i] = 3.0. From (2.2) it is found that any non-uniformity 
in concentration leads to a decrease in the observed viscosity, the magnitude of which is a strong 
function of the average concentration. Since the viscosity was observed to increase upon 
shearing, it was concluded that the initial concentration profile established upon loading the 
suspension into the gap was more non-uniform than that corresponding to Couette flow. 

2. 1 . Diffusion model for the observed short-term viscosity increase 

The viscosity increase phenomenon was modelled by assuming that the diffusion coefficient 
across the gap was constant throughout the migration. This approximation is acceptable if the 
variation in concentration is sufficiently small, and in any case the experiment yielded some 
average value of the diffusion coefficient for the migration. If the initial concentration profile is 
given by 0(0, y) and the concentration after a long period of shear is 0 (assumed to be constant), 
the concentration profile at all times is given by: 

0(y,t) = 0 + X B n cos ^ ex K“~T“0 (2,3) 

n= 1 b 

where 

Bn = C ° S ^” ^ 

and D,| is the diffusivity within the plane of shear parallel to gradients in fluid velocity. 

To further simplify the model, equation 2.3 was approximated by its limiting form at long 
times, i.e. the coefficients were chosen such that: 

B^fr, B n =0, n#l (2.4) 

in which case, to obtain accurate values of the diffusion coefficient under this assumption, 
equation 2.3 was fitted only to data taken after sufficient time has elapsed for the neglected terms 
to become unimportant. Since these higher-order terms decrease exponentially with a rate 
constant at least four times that of the leading-order term, this requirement was easily met Thus 
the model describing the short-term viscosity increase contains three adjustable parameters: the 
equilibrium uniform concentration; the amplitude of the initial variation in concentration Bj; and 
the diffusion coefficient. The first of these was fixed by the equilibrium viscosity corresponding 
to the concentration of the suspension initially loaded; thus the rate at which the viscosity 
approached its equilibrium value and the magnitude of the deviation between this and its value at 
the start of the experiment yielded the diffusion coefficient and the approximate initial 
concentration variation across the gap 

It is important to note that the above development tacitly assumes that the concentration 
profile is not a function of the distance up the gap. In practice this is unlikely to be true since at 
the base of the gap the concentration profile will correspond to the entrance region into the gap, 
while sufficiently far up the gap it should correspond to the steady-state distribution resulting 
from Poiseuille flow. While this will be discussed in more detail in section 4, it is noted here that 
while variations in the initial concentration profile up the gap may affect the estimated amplitude of 
the concentration variation Bj, it should not affect the calculated value of D N . This arises from the 
assumption that the diffusivity is essentially constant within the gap, corresponding to the average 
concentration which does not vary along the length of the gap, and thus the rate at which the 
concentration approaches its steady-state distribution is the same at all positions along the gap. 

2.2 Experimental results 

Short-term viscosity increase experiments were performed by Leighton and Acrivos (1987b) 
with suspensions of 46pm and 87|im polystyrene spheres at concentrations from 30% to 50% in 
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two mixture of silicone oils with viscosities of 1.22p and 1.07p. The measurements were taken 
using touette gap widths of 0.639mm, 1.261mm and 2.513mm. The gap height for all 
experiments was 4.508cm, and the bob diameter was 4.7498cm. 

The fit of the data to the model was excellent (cf. Figure 3); however, owing to the many 



Figure 3 Short-term viscosity increase with shearing, 46 pm polystyrene suspension. Model 
parameters were: <p = 0.45, y = 2.4 s 1 , — 1.8, B x — 0.062. 


assumptions that were necessary in deriving (2.3), the calculated values of the diffusion 
coefficient and concentration fluctuation across the gap must be considered only approximate. 

The diffusion coefficient was found to be proportional to ya 2 and, as in the case of diffusion 
normal to the plane of shear measured in the long-term viscosity decrease experiments, was a 
strong function of concentration. A plot of the diffusion coefficient as a function of concentration 
is given in Figure 4, where the dashed line is the diffusion coefficient normal to the plane of 



Figure 4 Diffusion coefficient observed in the short-term viscosity-increase experiments: O. 
46 urn polystyrene in medium gap; □, 46 pm polystyrene in small gap; A, 87 pm polystyrene in 
medium gap. Dashed line is the diffusion coefficient observed in the long-term viscosity -decrease 

experiments. 
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shear. The fact that the measured values of the two diffusion coefficients were almost identical 
provides us with added confidence that the viscosity increase was interpreted correctly as resulting 
from a particle migration across the gap width. 

The short-term viscosity increase phenomenon was observable only for a narrow range of 
particle diameters, suspension concentrations and Couette gap widths. This was due in part to the 
ract that, since the diffusion coefficient was found to be proportional to the square of the particle 
radius, the total length of shearing necessary to reach steady-state was inversely proportional to 
a , the square of the ration of the particle diameter to gap width. Thus, for very large values of ft 
(such as were obtained for experiments with the 87pm polystyrene spheres in the narrow Couette 
gap) the strain associated with the migration was sufficiently short that it was not possible to 
reliably separate fluctuations in the viscosity due to migrations across the gap from those due to 
theimtial equilibration of any short range order in the suspension, first observed by Gadala-Maria 
and Acnvos (1980). Similarly, the timescale for the viscosity increase was also much too short 
tor experiments with the 46pm spheres in the narrow gap and a 50% concentration, owing to the 
high value of the dimensionless diffusion coefficient found at this concentration. 

At 30% solids concentration, a different experimental difficulty was encountered in that 
although the total fluctuation in concentration across the gap may have been the same as that 
observed for more concentrated suspensions, the resultant variation in the observed viscosity was 
too small to be accurately measured. From the observed dependence of viscosity on 
concentration, the same fluctuation in concentration across the width of the gap affects the 
observed viscosity at an average concentration of 30% by an amount that is less by an order of 
magnitude than that at 50%. 3 

Finally, for some combinations of concentration, particle diameter and gap width the 
short-term viscosity increase effect was not observable. No useable measurements were obtained 

s '$* S . u |P ensi T, S of either 46 P m ° r S7 ^ m s P heres in ^ lar 8 e Couette gap at concentrations of 
30% to 45%, and for suspensions of 46pm spheres at 30% to 40% concentration in the medium 
Couette gap. Possible causes for the absence of a measurable initial increase will be discussed in 
section 4. Table 1 presents the calculated amplitude of the concentration fluctuation across the 
channel for those experiments where it was measurable. 


Particle 

Gap 




diameter 

width 




(pm) 

(mm) 

4> 

1 5 , 

B 1 

46 

0.639 

0.40 

l.i 

0.131 



0.40 

0.54 

0.094 



0.45 

2.45 

0.075 


1.261 

0.45 

1.7 

0.062 



0.50 

3.2 

0.081 

87 

1.261 

0.40 

1.4 

0.050 



0.40 

1.4 

0.075 



0.45 

1.6 

0.063 


Table 1 . Estimated diffusion coefficient and amplitude of concentration fluctuation for the 

short-term viscosity-increase experiments 
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3. Mechanisms leading to shear-induced migration 

Thus far, we have presented experimental evidence for the existence of shear-induced particle 
migrations arising from gradients in concentration and shear stress. To see how these migrations 
take place we follow the development provided by Leighton and Acrivos (1987b). Let us first 
examine the case of migration due to gradients in particle concentration. 

Consider a single marked sphere of radius a in a suspension of otherwise identical spheres 
undergoing the viscous linear shear flow u = yy, where u is the velocity in the x-direction. As the 
sphere interacts with its neighbors in the shear flow, it will experience a series of displacements in 
both the y- and z-directions with a characteristic length proportional to a and a frequency 
proportional to the shear rate y. In the absence of any gradient in concentration, these 
displacements will be random with zero mean (i.e. on average the particle will remain on its initial 
streamline), and thus will constitute a random walk. But, as is well known, in the presence of a 
concentration gradient such a random walk will lead to a diffusive flux and thus may be 
characterized by a diffusion coefficient, in this case with the dimensional scaling 7 a 2 , the same as 
was observed in the experiments described here. It is important to note that in a dilute suspension 
the spheres will return to their initial streamlines at the end of all two-particle interactions owing to 
the linearity of the viscous-flow equations when only viscous forces are present. As a 
consequence, at least three particles must interact to yield the permanent displacements that lead to 
a random walk, and therefore, since the reate at which two particles interact with the marked 
sphere is proportional to yj) 2 , the diffusion coefficient must be proportional to ffia 2 in the dilute 
limit. 

This source of diffusive flux, termed shear-induced self-diffusion, may be measured by 
examining the random walk of a single marked particle in a homogeneous suspension. 
Experiments carried out along these lines have been conducted by Eckstein, Bailey & Shapiro 
(1977), and more recently by Leighton & Acrivos (1987a), and have demonstrated that the 
diffusion coefficient is indeed proportional to ya 2 . Leighton & Acrivos (1987a) also obtained a 
dimensionless value of the diffusivity of about <J ) 2 /2 in the dilute limit, in agreement with the 
scaling predicted by the theory. The measured coefficient of self-diffusion, however, was found 
to be a much weaker function of concentration and, at a concentration of 40%, had a value nearly 
an order of magnitude lower than the diffusivity that was calculated from the experiments 
described earlier in this work. The discrepancy between the observations of self-diffusion and 
effective diffusivity in the presence of a gradient in concentration suggests, therefore, that the 
presence of a concentration gradient in some way induces a drift of particles from regions of high 
to low concentration in addition to that provided by random self-diffusion. 

The most likely source of this additional drift is that interparticle interactions in the presence 
of a gradient in concentration lead to an average displacement of the marked sphere from regions 
of high to low concentration. It is not clear at this stage whether interactions in the presence 
solely of viscous forces can lead to such a drift. Direct calculations of the drift for such a 
suspension would require consideration of interactions involving at least three spheres in a dilute 
suspension and, in concentrated suspensions where the average interparticle separation distance is 
very small, the interaction of many spheres would have to be taken into account. Such 
calculations are far beyond the capabilities of current analytical techniques, and are also quite 
difficult to deal with numerically owing to the very large number of particles that must be included 
in any computation. 

As we shall presently demonstrate, however, for sufficiently concentrated suspensions of 
real non-colloidal particles, it is not necessary or even appropriate to consider only the influence 
of purely viscous hydrodynamic forces because, in such systems, the particles are driven 
sufficiently close together by the flow that irreversible surface contact will occur as a consequence 
of surface-roughness effects, thereby destroying the macroscopic reversibility of the purely 
viscous interactions. In this section we shall therefore discuss both theoretical and experimental 
evidence for the existance of such irreversible interactions and demonstrate that they may lead to 
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the diffusivities observed here. 


3.1. Irreversible interactions in concentrated suspensions 

Consider a sphere interacting with a second sphere in a simple shear flow. As the two 
spheres approach one another, the viscous stresses in the fluid act to drive the spheres together. 
Under purely viscous conditions, when the interparticle separation distance is very small this 
approach is resisted by the lubrication layer between the particles, in which the resistance is 
inversely proportional to the separation distance. Thus although two mathematically smooth 
interacting spheres may never touch, hydrodynamic theory predicts that, over a certain range of 
initial configurations, they will approach one another very closely even in a dilute suspension. On 
the other hand, in the presence of a finite amount of surface roughness on the spheres (as is the 
case for any real particles), the interaction may be significantly modified. Indeed, Arp & Mason 
(1977) found that even for two isolated interacting spheres, the existence of a small degree of 
surface roughness was sufficient to eliminate the closed orbits predicted for purely viscous 
interactions. 

This effect of surface roughness is accentuated for concentrated suspensions. Specifically, 
since the forces driving spheres together in the flow depend on the bulk fluid stresses, they are 
proportional to the bulk suspension viscosity, which in turn is a strong function of concentration. 
In contrast, the lubrication forces resisting this approach remain proportional to the pure-fluid 
viscosity since the presence of other particles in the suspension does not affect the flow in the 
narrow gap between particles. This imbalance, combined with some finite surface roughness, 
implies that in sufficiently concentrated suspensions particles will simply approach one another 
without any significant displacemtn from the their original streamlines until they come into 
physical contact, Mowing which they will rotate owing to the vorticity of the shear flow, and 
finally separate. Moreover, since the interaction is no longer reversible owing to the surface 
contact it will also no longer be symmetric and thus will lead to permanent displacemtns of the 
particles from their original streamlines at the end of each interaction. A more complete 
discussion of the evidence for irreversible interactions and thier effect on the rheology of 
concentrated suspensions is given by Leighton (1985). 

3.2 Particle drift arising from irreversible interactions 

There are several ways in which the irreversible interactions described above can lead to drift 
in the presence of gradients in concentration or shear stress. To see this, consider a test sphere at 
the origin which is immersed in a suspension undergoing shear. We shall assume that the bulk 
flow is in the x-direction with a constant shear stress a = a yx and a linear concentration gradient 
in the z-direction, normal to the plane of shear. Under these conditions, the viscosity, and hence 
the shear rate, will be be constant within the plane of shear. Consequently, when the particles are 
not mathematically smooth, the test sphere will be irreversibly displaced upwards following an 
interaction with another sphere approaching it from below, and conversely if approached from 
above. Thus, in the presence of a higher particle concentraiton on one side of the test sphere than 
on the other, this sphere will experience a drift towards the region of lower concentration since it 
interacts with more particles on one side than on the other. Moreover, the displacement after each 
interaction will be proportional to the particle radius; thus since the excess rate of interactions from 
regions of higher concentration is proportional to vad<|)/dz, the particle flux resulting from this 
source of drift should be proportional to 'ya 2 <{xi(J)/dz, which is the same scaling expected for 
shear-induced diffusion. 

A second source of drift normal to the plane of shear arises from the gradients 
in suspension viscosity brought about by gradients in concentration. First we note that in the 
absence of a gradient in viscosity, two touching spheres in a shear field will rotate about their 
center of mass (the point of contact). But, in the presence of a viscosity gradient, the center of 
mass will no longer be the center of rotation and the particles will, on average, be displaced 
during an interaction from regions of high to low viscosity. The magnitude of this displacement 
during each irreversible interaction will scale as the relative variation in viscosity across the 
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particle, i.e. it will be proportional to (a/|i)d|i/dz, multiplied by the particle radius. Since, for 
small gradients in concentration, the variation in viscosity is linear in the concentration gradient, 
the above result, when multiplied by the rate of interactions ytj>, gives the corresponding drift 
velocity. In turn, when multiplied by 4> and divided by d<|>/dz, the drift velocity leads to an 
expression for the effective diffusivity arising from this source of drift D ~ ya z (<^/|i)d|Vd<]>. The 
total effective diffusivity normal to the plane of shear is then the sum of the two mechanisms 
outlined above, plus that due to shear induced self-diffusion. 

For concentrated suspensions, the function (l/|i)d|id<{> is very large, with a value of about 15 
for a 45% suspension; thus drift due to gradients in viscosity is likely to dominate the diffusivity. 
Recognizing this, we therefore let 


n v * d ^2 

Di " Kl F5 1 “ 


(3.1) 


be the expected form for the effective diffusivity normal to the plane of shear at high 
concentrations, where K_l is a dimensionless parameter whose value will depend on the exact 
geometry of the interactions. As shown in figure 5, Kj_, as determined from the results of the 
long-term viscosity experiment, was found to be a relatively weak function of concentration with 
a value of about 0.7 at 45%. The scatter in the data is, of course, indicative of the simplcity of the 


The contribution of irreversible interactions to drift within the plane of shear is quite similar 
to that found above for the case of drift normal to the plane of shear. Again, displacements will 
lead to random self-diffusion, to drift arising from a higher rate of interactions on one side than on 
another due to the concentration gradient (assuming a uniform shear rate), and drift from regioins 
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Figure 5 Value of K L vs. concentration: □, 46 pm polystyrene in 0.639 mm gap; O, 46 pm 
polystyrene in 1.261 mm gap; A, 87 pm polystyrene in 1.261 mm gap. Open symbols, y - 24 s ; 
filled symbols, y = 76 s' 1 ; half-filled symbols, y = 7.6 s l . 


of high to low viscosity. In addition to these sources of drift, however, variations in viscosity 
within the plane of shear will, in the case of a uniform applied shear stress, lead to variations in 
the local shear rate. Thus, in regions of low concentration and low viscosity the shear rate will be 
greater, with the consequence that a test sphere will, on average, experience a greater number ot 
interactions from the region of lower concentration than would otherwise be the case. Since the 
shear rate is inversely proportional to the local viscosity for a constant shear stress, the excess rate 
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of such interactions will be proportional to (<Jrya/|j.)(dp/dy), which in turn will reduce the drift 
Y. ® re S lons high to low concentration by an amount proportional to 

Cya ^(dp/d^jCd^/dy), an expression for the drift velocity identical with that due to gradients in 
viscosity with uniform shear. 

, AUhough this drift due to the shear rate gradient effect will certainly reduce the magnitude of 
the effective diffusivity (the sum of all contributions due to drift and random walk) it appears 
unlikely that me two terms due to gradients in fluid viscosity will exactly cancel out. Thus as in 
the case of diffusion normal to the plane of shear, we obtain for the limiting expression of the 
effective diffusivity at high concentrations 


n _ ^ dp . 2 

Di= Ki —ya z 


H d<j> 


( 3 . 2 ) 
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where the value of K ( | may be determined from the initial- viscosity-increase experiments. The 
results are shown in Figure 6 from where it is seen that, in spite of the scatter, which again is 
indicative of the simplicity of the model, K„ appears to be relatively independent of concentration 
and approximately equal to 0.6. 

• . • Th ? sai ? e mech anism due to gradients in shear rate that reduced the effective diffusivity 
within the plane of shear should also lead to drift from regions of high to low shear stress in a 
homogeneous suspension. This is because, for the case of uniform particle concentrations and 
uniform suspension viscosities, the local shear rate is proportional to the local shear stress <r 
hence the excess rate of particles ; interacting with the test sphere from regions of high shear stress 
is proportional to (<jrya/G)/(dG/dy), yielding a particle flux 

2 

<t> da 


Ny=-K C - 


dy 


7a 


( 3 . 3 ) 


where again K 0 is some order-one function of concentration. 



Figure 6 Average value of K, vs. concentration calculated from the data in table 1 
denote one standard deviation as estimated from the scatter in the experimental 


• Error bars 
results. 
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4. Particle distributions in Poiseuille How 

Using the above analysis it is possible to determine the steady-state concern ration ^mbunon 

gsssiissss?- 

particle flux: 2 
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da „ i_dM.^1V 
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o ay " d<J) dy- 

is set equal to zero. For Poiseuille flow the shear stress is simply proportional to y, hence 

I do _ 1 (4.2) 

a dy y 

and therefore we obtain the derivative of the steady-state concentration pro e. 

^ = 3a /idMV 1 ! (4.3) 
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profile: 
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(4.4) 


iL= ^K, 

may uivert equation 4.4 to obtain the concentration profile across the channel. 
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K^j/Kii is g rj Tn,:„ a conseauence of the singularity in (4.2), and the 

S^reonhe^effKhve diffasivity as the particle concentration approaches its maximum value. 

4 ‘‘ wemKsfa* the value of from the amplitude of the short-term viscosity increase 

phenomenon if we examine in more detail the procedure by winch wn»»n 

rtpvirp As is described in detail by Leighton and Acnvos (1987b), in these experiments 
the fluid was loaded into the gap by first pouring the suspension into the Couette cup and then 
SwerinX boMhuTSspfS the fluid and filling the gap. The flow thus created cons, sts of a 
converging entrance flow at the base of the gap, followed by channel flow up die gap. Wide it « 
nnelear what effect if any, the entrance flow has on the concentration profile, the channel flow up 
the eao should lead to a migration of the particles from the regions near the walls into the center. 
Theconcentration profile in the gap will thus be a function of distance up the gap, but provide 
22 SS heff of the gap is sufficiently large the concentration profile in the gap will approach the 

steady-state distribution given by (4.5). 
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Figure 7. Plot of concentration profile across a channel predicted from equation 4 5 for an average 

concentration of 45%, as a function of K C/ / K , . S 
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resulting in the dimensionless equation: 
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which may be solved using separation of variables. Equation 4.8 admits a solution 


of the form: 
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(4.9) 



While equation 4.10 does not admit a closed form solution, the leading eigenvalue has been 
determined numerically to be ocq = 2.277795. As a consequence, the assumption that the 
concentration distribution in the gap is close to that at steady-state for channel flow will be valid 
provided 

2 <xg a^6,K >> , (4 . n) 

where h is the gap height. In the experiments described by Leighton and Acrivos (1987b) where 
the short-term viscosity increase was observable, the dimensionless parameter given in (4.1 1) 
ranged from 1.9 to 14, thus the assumption of steady-state is reasonably good. In contrast, in 
those experiments where the short-term viscosity increase was not observed, the value of the 
parameter was less than 1, suggesting that the observed increase was, in fact, due to particle 
migrations in Poiseuille flow rather than entrance effects. 

To estimate the value of K„, we may simply substitute the steady-state concentration 
distribution given by equation 4.5 into the equation for Bj given in equation 2.3 and integrate. 
The resulting values of K 0 are given in Figure 8, which, allowing for considerable scatter in the 
experimental data, indicates that is a relatively weak function of concentration (not quite 
constant as we have assumed) with a value of about 0.6 at <t> = 0.45. 
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Figure 8 Average value of K a vs. concentration calculated from the data in table 1. Error bars 
are one-standard-deviation error estimated from the scatter in the experimental results. 
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4.2 Comparison with other experiments 

In a bounded Poiseuille flow (channel or tube), any migration of particles toward the center 
leads to an increase in the viscosity in this region and hence to a blunting of the parabolic velocity 
profile that applies for Newtonian fluids. We may calculate the expected velocity profile resulting 
from the steady- state viscosity distribution given in equation 4.4: 



which, surprisingly, is not a function of concentration. Of course, in the development leading up 
to (4.4) we assumed that the ratio K^ZR,, is independent of concentration. This is certainly not 
true for dilute suspensions where the mechanisms leading to the shear stress gradient induced 
migrations would be expected to vanish, and is only approximately correct for concentrated 
suspensions. 

Observations of velocity profiles in concentrated suspensions flowing through tubes were 
conducted by Kamis, Goldsmith and Mason (1966). Their observations at an average 
concentration of 38% (the highest concentration for which measurements were reported) and 
particle/tube diameter ratio of .028 is reproduced in Figure 9 together with the profile predicted by 



Figure 9. Comparison of the velocity distribution estimated from equation 4. 12 with that observed 

by Kamis, et al. (1966). 


equation 4.12 using an observed ratio B^/K,, = .912 estimated from the 0 = 0.40 experiments 
described above. From this comparison it is seen that the velocity profile observed by Kamis, et 
al. is blunter than that predicted here, corresponding to Kq/Kh = 2.1. While these two values do 
not greatly differ considering the large degree of scatter in our experiments, at least part of the 
discrepancy may be accounted for by wall effects. Indeed, Kamis, Goldsmith, and Mason 
attributed their observed blunting entirely to wall effects, however in view of the small particle 
diameter / tube diameter used in the experiment depicted here, this seems unlikely. Experiments 
in the Couette geometry reported by Kamis, et al. (1966), showed that wall effects were 
insignificant at the same concentrations even for particle diameter / gap width ratios over a factor 
of three larger than that used in the experiments leading to Figure 9. 
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These authors did attempt to measure the concentration profile across a tube, although only at 
a lower average concentration of 32%. Owing to the statistical nature of their measurement 
technique, however, their observations had a very large experimental error. Specifically , their 
measurements involved counting the number of spheres that were present in each of four 
divisions of the half- width of a tube cross-section, and since only a small fraction of the spheres 
were marked, their observations were subject to Poisson statistics, which dictate that the 
one-standard-deviation error in the number of spheres that were counted is equal to the square 
root of that number. But since the total number of spheres that were counted in each region was 
rather small (less than 80), the concentrations reported by Kamis et al. (1966) had a two standard 
deviation error of nearly 25%. Their observations, with error estimated as described above, 
together with the concentration profile predicted using equation 4.5 and a value of Kg/Kn = .912 
is given in Figure 10. While the concentration observations of Kamis et al. do not agree with the 
predicted concentration profile (which has been extrapolated well beyond the range over which 
Ko/Km was estimated), more accurate measurements of the concentration distribution in Poiseuille 
flow are needed in order to determine the source of the observed blunting of the velocity profiles. 



Figure 10. Comparison of concentration distribution estimated from equation 4.5 with particle 
distribution observed by Kamis, et al. (1966). Error given is 2 standard deviation error 

calculated by statistical means. 
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Abstract 


Liquid jet breakup mechanisms and processes 
have been studied extensively over the last one 
hundred years. However, since the region near the 
jet injector is too dense and optically opaque, 
conventional visualization cannot be applied with 
satisfactory experimental results. To unravel the 
liquid Jet breakup process in the non-dilute 
region, a newly developed system of real-time 
X-ray radiography together with an advanced 
digital image processor and a high-speed video 
camera was used in this study. Based upon 
recorded X-ray images, the inner structure of a 
liquid jet during breakup wa 3 observed. The Jet 
divergent angle, jet breakup length, and void 
fraction distributions along the axial and 
transverse directions of liquid jets, etc., were 
determined in the near-injector region. Both 
wall- and free-jet tests were conducted to 
study the effect of wall friction on the jet 
breakup process. 


Introduction 


Understanding the breakup mechanism of liquid 
jets in diesel engines, regenerative liquid 
propellant guns, and many other combustion and 
propulsion systems is extremely important, since 
the jet breakup process can strongly affect the 
jet divergent angle, the distributions of droplet 
size and velocity, and the mixing and combustion 
processes. As indicated in a detailed review 
paper by Arcoumanis and Whitelow, 1 numerous 
theoretical and experimental attempts have been 
made to study the stability and breakup phenomena 
of liquid jets ejected from nozzles. In recent 
years, Wu, Reitz and Bracco, 2 Birk and Reeves, 3 
and Baev et al.^ conducted intensive studies in 
the experimental investigation of breakup behavior 
in the near-injector region. Shimizu et al. 5 
studied the breakup length of a high-speed Jet by 
measuring the electrical resistance between a 
nozzle and a fine wire detector located in a spray 
jet. They found that breakup length decreases 
with increase in injection velocity, finally 
reaching a constant value. 

In the theoretical analyses, the growth of 
initial pertur bations on the liquid surface, the 
effects of liquid inertia, surface tension, 
viscous and aerodynamic forces on the jet were 
considered. The earliest development of a 
predictive model for the jet breakup was initiated 
by Rayleigh. ^»7 His linear stability analysis of 
an inviscid cylindrical liquid jet showed that an 
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axisymmetric disturbance could be stable or 
unstable, depending upon the magnitude of the 
wavelength (stable for wavelengths less than the 
circumf erence of the jet, and unstable for other 
cases). Further investigations using linear 
stability analysis were conducted by Tyler and 
Richardson, ^ Schweitzer. 9 Merrington, 1 ^ Levich 1 1 , 
Dombrowski and Hooper, 12 and others. Based on 
these studies, the Jet divergent angle near the 
injector exit, and jet breakup length and mean 
drop size near the lateral surface of the jet can 
be estimated from the Reynolds number, the Weber 
number, the density ratio of liquid to gas, and 
several empirical constants. The effect of 
viscosity was studied by Weber . 1 3 More recently, 
the effect of nonlinearity 1 ^ -1 9 has been 
considered and solved numerically in order to 
simulate more closely the jet breakup processes. 

In the experimental studies, the core of the 
jet is optically opaque due to the high-density 
condition in the near-injector region. Hence, 
observations and measurements of the near-injector 
region are almost impossible. Consequently, most 
previous studies focused on the dilute region or 
on the jet boundary near the injector exit. 2 Due 
to the lack of quantitative test results in the 
near-injector region, the theoretical models 
developed for prediction of jet-core breakup 
length, mean drop size distribution, etc., have 
not yet been validated. The estimated jet 
profiles and characteristics at a given 
near-injector station were used as initial and/or 
boundary conditions in order to simulate the 
atomization and combustion processes in the dilute 
region; this, in turn, could cause errors and 
uncertainties in the theoretical predictions of 
combustion processes. 

The purpose of the present study is to 
achieve a better understanding of the jet breakup 
process, using a newly developed system of 
real-time X-ray radiography together with an 
advanced digital image processor and a high-speed 
video camera. The specific objectives are: (1) 

to illustrate the differences between X-ray 
radiography images and regular high-speed movie 
films; (2) to measure jet velocity, jet breakup 
length, and void fraction distributions in the 
near-injector region; (3) to observe the evolution 
of the jet breakup and the formation of ligaments 
and droplets in the near-injector region in order 
to determine the inner structure of the liquid 
jets; and (4) to demonstrate the feasibility and 
advantage of using X-ray radiography in the jet 
breakup study. In order to study the effect of 
wall friction on the Jet breakup process, both 
wall- and free-jet tests were conducted. 

Experiment Apparatus 


In the experimental approach, a test rig has 
been designed and fabricated to simulate both wall 
and free jets ejected from a two-dimensional slit 
with an aspect ratio of 10. The advantage of 
using a planar two-dimensional jet is to avoid the 
curvature effect which is inherently associated 
with circular or annular jets. 3 The side view of 
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the jet breakup process is much clearer in the 
planar two-dimensional case than in other cases. 
The schematic diagram of the test rig is shown in 
Fig. 1 . A photograph of the test setup is shown 
in Fig. 2. In Fig. 1 , the gap of the 
two-dimensional slit between elements 5 and 7 is 2 
mm, and element 5 is an interchangeable piece for 
providing different lengths of the wall jet. When 
a piece of element 5 ends at the same location as 
element 7 , the configuration is that of a free jet 
test case. A pressure transmitting rod is shown 
as element Prior to the test, this rod was 
lifted by a vacuum pump. The liquid was then 
loaded through element 6 into the free volume 
beneath the rod and the entire slit region. The 
liquids used in tests were either water or 
Pantopaque (C 19 H 29 IO 2 )* a fl uid used in medical 
X-ray examinations. The physical properties of 
the two test fluids are presented in Table 1. 
During the test, a solenoid valve was activated to 
introduce high-pressure nitrogen gas into the top 
portion of the test rig (see Fig. 2), thus pushing 
the rod downward to inject the liquid into the 
unconfined ambient air. A major portion of the 
Pantopaque was recovered for further use. 

SOLENOID VALVE 



Fig. 1 Schematic Diagram of Test Rig Assembly 
for Two-Dimensional Wall or Free Jet 
Breakup Studies (1. top cover, 2. bolts, 
3 . liquid reservoir, 4. push 
rod, 5 . top plate of slit, 6 . liquid 
feeding hole, 7 . base plate of slit, 8 
and 9 . test rig holder) 



Fig. 2 Photograph of the Test Setup 


Data Acquisition System 

During the test, the instantaneous liquid jet 
contour was filmed by real-time X-ray radiography. 
Figure 3 shows the layout of various components of 
the radiography system. A continuous X-ray was 
generated from the Phillips MG 321 constant 
potential X-ray system. Two X-ray tubes [one with 
a focal spot combination of 0.2 x 0.2 mm/ 3.0 x 3,0 
mm (MCN 167/160 kV) , the other with a focal spot 
combination of 1.2 x 1.2 mm/4.0 x 4.0 mm (MCN 
321/320 kV)] were employed to achieve different 
penetration depths and spatial resolution 
requirements. A lead diaphragm was installed at 
the exit port of the X-ray tube head to limit the 
angle of divergence of the X-ray beam and to 
confine the beam to the measuring section of the 
liquid jet test rig. (This also reduced 
unnecessary radiation exposure.) A second lead 
diaphragm with a larger opening was placed in 
front of the image intensifier to reduce scattered 
X-ray radiation and decrease the noise level on 
the fluorescent screen. 

After passing through the test rig and liquid 
jet, X-ray signals were transformed to fluorescent 
light signals on the output screen of a tri-field 
image intensifier (Precise Optics, Model PI 2400 
ATF, 4 " , 6 " , or 9" field diameter). The input 
fluor of the image intensifier was made of cesium 
iodide with a decay-time constant of 650 ns, and 
the output fluor was a p 20 type with a 85 ns 
decay-time constant. These time constants are 
short enough to allow the motion analysis system 
to operate at its maximum framing rate without 
generating image blur. 

The fluorescent light signal output from the 
image intensifier was recorded by a Spin Physics 
2000 Motion Analysis System. This system can 
record up to 2000 fully digitized frames per 
second, or up to 12,000 digitized pictures per 
second with adjustable playback speed. The motion 
analysis system consists of the following 
subsystems: 

a) a Spin Physics 2000 video camera with 
solid-state image sensor. The picture 
information goes to the console from 
the camera, and is processed into a 
frequency-modulated carrier that is 
recorded on tape. 

b) a main electronic bin, which contains 
record and playback electronics, along 
with video output circuitry. 

c) a tape transport, which drives the 
half- inch tape cassette at a maximum 
speed of 250 inches per second. 

Digitized data was stored on a high-intensity 
magnetic recording tape and transferred to the 
digital image-processing system frame by frame for 
analysis through an IEEE-488 interface. The 
digital image-processing system consists of 
several major components: 

a) a Quantex (QX-9210) digital image 
processor with two pipeline point 
processors, each with a random access 
image memory of 480 x 640 x 12 bits. 

The processors perform real-time image 
enhancements, including noise 
reduction, image subtraction, arbitrary 
contrast control, roam, and zoom. The 
processor can also be used to conduct 
such high-speed image analyses as 
brightness histogram, local contrast 
stretch, area and point brightness 
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Fig. 3 Real Time X-ray Radiography and Digital 
Image Processing System 

measurements, calibrated length and 
area measurements, Sobel edge 
enhancement, and so forth. 

b) a mass storage device consisting of 
both a 2.4 MB dual 8" floppy diskette 
drive and a 50 MB hard disk for storage 
of image data. 

c) a high-resolution videograph printer 
(AIDI CT1500) for producing large 
high-quality pictures on hard copies 
with 1660 lines/inch resolution. 



X 


X 


Fig. 4 X-ray Intensity Distribution Across a 
Liquid Jet Under Idealized Conditions 



Data Reduction Procedure 

Besides observing the Jet structure directly 
from the recorded film (both regular movie and 
X-ray films), more detailed and accurate data 
could be deduced by analyzing the X-ray intensity 
distribution across the Jet in horizontal and 
vertical directions. Basically, an ideal 
radiography image of the liquid Jet can be 
determined from the assumptions that 1 ) X-ray 
radiation is generated from a point source via an 
infinitely small focal spot, 2) X-ray is only 
attenuated by photoelectric absorption, and 3) 
distribution of X-ray intensity over the input 
screen of an image intensifier is uniform. As 
shown in Fig. 4, the X-ray intensity distribution 
on a y - constant plane can be evaluated by the 
following equation: 

I(x) - 1 0 EXP [-mW(x) ] (1) 

where m is the linear attenuation coefficient of 
photoelectric absorption, and W(x) is the local 
sum of the intercepted width of the liquid in the 
two-phase jet. Due to the X-ray beam attenuation 
across the liquid jet, the following test data can 
be obtained using the digital image-processing 
system: 

a) Measurement of Radiance or Pixel Value: 
The radiance (or pixel value) of any 
arbitrary point of interest in the jet 
can be obtained by moving the cursor to 
that point. The local radiance and the 
x, y coordinates of the location can be 
displayed on the screen. 

b) Measurement of Area and Pixel Value: 

The software program combines a 
measurement of an area with a 
measurement of the integrated pixel 
values (or total radiance) within the 
area. The calculated average pixel 
value within the area is also displayed 
on the screen. This feature can be used 
to detect the local void fraction level 


from the X-ray radiograph for a 
specified jet region (see Fig. 5). 

c) Profile Analysis: Figures 6 and 7 show 

the measurement of the pixel value of 
two scan lines — one parallel and one 
perpendicular to the jet core and head 
regions. Both end at the cursor point. 
From these profiles, the void-fraction 
distribution can be deduced for a 
certain cross-section of the jet. 

d) Histogram Analysis: This procedure 

consists of measuring the frequency of 
occurrence of pixel values within a 
selected boundary of a rectangle. Such 
distributions can be used to detect jet 
breakup phenomena in the core region. 
Keeping the axial length of rectangular 
area constant by varying the width of 
the rectangle, different sets of 
histograms can be obtained. These are 
arranged into a series of histograms 
(see Fig. 8). 

e) Isophote Analysis: Using the data 

reduction program, certain regions of 
the viewing area with the same range of 
pixel values can be replaced by white 
spots. This feature makes it possible 
to observe the contours of different 
radiances. The inner structure of a 
liquid jet can, therefore, be studied by 
selecting different pixel values (see 
Fig. 13, to be discussed later in the 
results section). 

Finally, jet-head velocities in both axial 
and radial directions can be deduced from the 
recorded film by using the elapsed time and 
displacement of reticle lines on the monitor of 
the Spin Physics Camera System. 

Discussion of Results 

A series of high-speed motion pictures 
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Fig. 5 Measurement of Area and Pixel Value of a 
Liquid Jet 



Fig. 6 Horizontal X-ray Intensity Profile at the 
Center Plane of a Liquid Jet 



Fig. 7 Vertical X-ray Intensity Profile at a 
Distance from Jet Exit 


showing the breakup process of a plane wall jet is 
given in Fig. 9. Several interesting phenomena of 
wall jet evolution and breakup can be noted from 
these photographs. In Fig. 9a, as the liquid jet 
traveled along the wall surface in the early 
stage, the head region diverged in the upward 
direction. A thin sheet of liquid was formed 
ahead of the main jet, as shown in Fig. 9b. This 
thin sheet of liquid surface is believed to be 
caused by shedding phenomena introduced by the 
inertia of the fast-moving liquid near the free 
surface of the engulfing wave front. Near the top 
surface, the velocity of the liquid is much higher 
than that near the wall, since the viscous force 
i3 much smaller at the free surface. As time 
progressed, this sheet became more evident, and 
the head region expanded further (see Fig. 9c). 

In Fig. 9d, the main jet accelerated and merged 
with part of the thin liquid sheet. In the later 
stage, the so-called Klystron effect 20 was visible. 
This effect was introduced by the acceleration of 
the main jet, which had higher axial momentum to 
catch up with the precursor jet mass. Due to the 
coalescence of the fast and slow moving fluids, 
the width of the jet spread in transverse 
directions. This spreading also made the jet-head 
portion more symmetric. In the meantime, several 
ligaments and numerous droplets were formed as the 
head region expanded further. Such a breakup 
process greatly influenced vaporization, ignition, 
and combustion of liquid sprays. 

A set of high-speed photographs for the free 
jet case is shown in Fig. 10. The jet head was 
quite symmetric as it expanded in the transverse 
direction. The shape of the jet head changed from 
a round cross section to a mushroom-shaped contour. 
The Klystron effect can also be seen in Figs. lOd, 
e, and f. The jet head in Fig. 1 0e became spear 
shaped with a sharp leading edge. Some ligaments 
and droplets were formed (see Fig. lOf) near the 
boundary of the liquid-gas interface. 

To demonstrate the differences between 
conventional high-speed movie film and high-speed 
X-ray radiography, a set of X-ray movie films 
showing free jet breakup processes is given in 
Fig. 11. It is quite obvious that, unlike the 
conventional high-speed movie film, the radiance 
of the jet head region is highly nonuniform. The 
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Fig. 8 


Histogram Analysis on 
Area of a Liquid Jet 


Different Selected 
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Fig. 9 Interesting Phenomena of a Plane Wail 

Liquid Jet Breakup Process (From 10 Evolution of a Plane Free Liquid Jet 

High-Speed Movie Film) Breakup (From High-Speed Movie Film) 


inner structure of the jet can, therefore, be 
observed and analyzed from these images. 

An enlarged photograph of the X-ray 
radiography, obtained from the zoom feature, is 
shown in Fig. 12. Figure 13 gives a set of 
pictures of the same object, but with different 
levels of radiance. The variation in contour of 
the isophote shows the detailed structure of the 
nondilute liquid jet. Although some pulsed X-ray 
photographs were obtained by Baev, et al., 4 for 
liquid jet, the detailed inner structure of a 
non-dilute jet during the process of breakup has 
never been observed before, to the best of the 
authors 1 knowledge. 

Using the Quantex Image Analyzer, the void 
fraction distribution in the vertical direction 
(normal to the jet axis) can be deduced for any 
axial location. Figure 14 shows the void fraction 
distribution across the widest jet head of 
photograph 12 (see also Fig. 7 cursor station). 
Near the centerplane (y * 0), the void fraction is 
much lower than that at the jet boundary [when y 
is approximately equal to 4H (gap width of the jet 
exit)]. The magnitude of the local void fraction 
is taken to be proportional to the local pixel 
value, (see Fig. 7 ) i.e., 

<J>(x,y) -It Ln 1 (2) 

where I G and I L represent the intensity (pixel 
value) of pure gas and liquid, respectively. For 
Pantopaque in air, the values of Iq and 1^ were 
130 and 20, corresponding to the X-ray setting ‘'or 


most tests. It is useful to note that the 
distance between the X-ray tube head and the image 
intensifier was held constant during these tests. 

A typical void fraction distribution for a 
wall jet is shown in Fig. 15. It is evident that 
near the wall (y/H » 0.5) the void fraction was 
almost zero, while the value of <p was quite close 
to unity near the outer edge of the wall jet. 

A question about the effect of diffraction 
and scattering from multiple interfaces on the 
accuracy of Eq. (2) has been addressed. Equation 
(2) essentially assumes that diffraction and 
scattering from multiple interfaces have no effect 
on the reduction of void fraction from intensity 
measurements. To verify the weak dependence of 
intensity on the number and orientation of 
interfaces, two strands of rectangular-shaped 
solid propellant grains from the same batch with 
equal length were tested. One strand was cut into 
10 sections with different angles to the main axis. 
The other strand remained as an integral piece. 

The two strands were placed on a platform adjacent 
to each other and viewed by X-ray from their end 
surfaces. The local intensity distribution for 
each strand was obtained and compared. It was 
found that the intensity of the uncut piece is 
only 5.80? higher than that of the chopped strand. 
This implies that the intensity of the X-ray 
radiography is minimally affected by the number 
and orientation of the interfaces. The loaded 
fraction (l-<t>) deduced 'rom Eq. (2) corresponds to 
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Fig. 11 Evolution of a Plane Free Liquid Jet 
Using the X-ray Radiography System 

the total width of the liquid intercepted by the 
X-ray. 

The jet spreading angle in the y direction is 
of interest in spray combustion. The angle 9 can 
be determined from the instantaneous Jet 
velocities in both x and y directions. The jet 
head velocities (u and v) in the x and y 
directions are deduced from the displacements of 
the vertical and horizontal lines on the Spin 
Physics monitor. The spread angle is defined as 
6 - 2 tan " 1 (v/u) (3) 

A plot of half spread angle (0/2) versus jet head 


Fig. 13 Isophote Analysis Shows the Inner 

Structure of a Liquid Jet Using Different 
X-ray Intensity Levels 

axial velocity (u) is shown in Fig. 16. For both 
wall- and free-Jet experiments, spread angles 
decrease as jet head velocity increases. This is 
due to the fact that when a jet has higher axial 
momentum, the rate of spreading is lower. Based 
upon the data shown in Fig. 16, the wall jet 
spread out more than the free jet due to the 
momentum distribution effect in the wall jet. 

Near the wall surface, the non-slip condition must 
be satisfied, and hence the liquid near the free 
surface has a higher velocity than the free jet 
case for the same Jet head velocity. 



Fig. 12 Enlarged Picture of a Liquid Jet from 
X-ray Image 
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The jet surface breakup distance (L SB ) and 
the jet core breakup distance (Leg) from the exit 
station are also important in spray combustion. 
Figure 17 plots these two distances versus jet 
head velocity. Determination of L BB was based 
upon the location of the first discernible 
divergence of liquid from the jet surface. The 
value of Leg was determined from the abrupt 
increase of intensity along the axial intensity 
distribution at the centerplane of the Jet. For 
both wall and free Jets, L sg increased with jet 
head velocity due to the higher inertia of the 
ILquid jet which delays the jet surface breakup. 
However, Lq B decreased while jet head velocity 
increased for both wall and free jets due to the 
fact that higher velocity jets may cause droplet 
collision and coalescence to occur earlier, 
generating an inhomogeneous axial intensity 
profile. It should be noticed that the jet core 
breakup distance used here is different from that 
used in previous studies. (In earlier studies, 

L C g is the length beyond which the Jet is no 
longer continuous.) At the same Jet head 
velocity, wall jets have longer L SB and L CB than 
free jets due to the fact that wall jets, which 
are affected significantly by the stronger viscous 
force, retard the core breakup. In comparing 
surface breakup length, at the same value of jet 
head velocity, the wall Jet has lower momentum 
than the free jet; hence, breakup length is 
longer . 

Summary and Conclusions 

1 ) The advantages of using real-time X-ray 
radiography for liquid- jet breakup 
measurements are summarized below: 

a) X-ray can partially penetrate through the 
liquid jets, even in the dense regions. 
Therefore, it is feasible to determine the 
instantaneous inner structure of the jet 
in the near-injector region during the jet 
breakup processes. 

b) X-ray radiography is a nonintrusive 
technique which does not affect the liquid 
jet breakup processes and is superior to 
other methods. 

c) The real-time feature of X-ray radiography 
gives the entire history of the jet 
breakup event instead of a few snapshots. 
The X-ray motion pictures taken during the 
test event can be played back immediately 
at a lower speed for detailed flow 
visualization and analysis* 

d) All data are in digital form and are 
convenient for recording, analysis, 
transfer, and storage in computers. 

2) Based upon the recorded X-ray images analyzed 
on a digital image processor, jet divergent 
angle, jet breakup length, and jet void 
fraction distributions can even be measured in 
the near-injector region. 

3) The effect of wall friction on the jet 
breakup process was observed by conducting 
wall- and free-jet tests. While the 
spread angle is larger for the same jet 
head velocity, longer distances are 

required for the surface and core to break 
up in the case of wall jet. 

4) To understand the jet breakup mechanism, 
more detailed studies are required. 

Future tests should be conducted under 
varied conditions, such as different 


pressure and density levels in the gas 
phase, higher liquid-jet velocity ranges, 
various test liquids, different jet nozzle 
geometries, etc. The development of a 
comprehensive theoretical model for the 
near-injector region is also necessary. 

The predicted results should be compared 
with experimental data for model 
validation. 
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TABLE 1 . PHYSICAL PROPERTIES OF TEST LIQUIDS 


IZE£ 

Water 

Pantopaque 

(C i 9H29IO2) 


Viscosity 

(poise) 

0.01 

0.051 


Density 

(gm/cm3) 

1.0 

1.259 


Surface 

Tension 

(dynes/cm) 

72.7 

32.1 



Fig. 14 Void Fraction Distribution Across the 
Widest Jet-Head Region (Plane Free Jet 
Case) 



Fig. 15 Void Fraction Distribution Across the 

Wildest Jet‘Head Region (Plane Wail Jet 
Case) 
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Fig. 16 Liquid Jet Divergent Angle Under 
Different Jet-Head Velocities 
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Editor’s Narrative: 


Round-Table Discussion: Physical, Model, Numerical and Experimental Aspects 
of Mixing and Demixing Processes in Multiphase Flows 


At the conclusion of the first day of presentations a discussion was opened 
in the general areas of physical, model, numerical and experimental aspects of 
mixing and demixing processes in multiphase flows. Two specific topics became 
the focus of the discussion. These were the relative merits of the continuum vs. 
Lagrangian-Eulerian frameworks for developing multiphase flow theories and the 
accuracy of the turbulence models used within these theories. 

By the lack of discussion, it would appear that laminar flow about a single 
spherical particle, at small or zero Reynolds number, leading to expressions for 
the drag and migrations in idealized mechanical two-phase flows are not imperative 
issues. This is despite the fact that most theories of turbulent multiphase flow, 
whether continuum or Lagrangain-Eulerian in origin, accept aspects of Stokes 
flow analysis within the models for the phase coupling drag. Hence, the accurate 
description of the forces on a single particle in various circumstances is important 
for a wide variety of flows. In addition, there exists a prevailing attitude that 
the diffusion potentials due to the turbulence of the continuous fluid phase 
dominates any transverse or migrational motions of the dispersed phase. How- 
ever, it should be noted that laminar processes resulting in transverse particle 
migrations may be of consequence within the viscous sublayer of turbulent flows 
if the particle sizes are small relative to the viscous sublayer thickness. This may 
be the case for such phenomenon as near wall particle depletion, as observed in 
turbulent multiphase pipe and channel flows. Lastly, within this same context, it 
should be noted that there are particle migration processes which may be ana- 
lyzed by in-the-mean or laminar theory for some scale of ’large’ or ’heavy’ 
particles relative to the turbulence. For example, a baseball, by imparting a spin 
to it can be made to curve (migrate laterally) through the atmospheric boundary 
layer with apparent disregard for the random diffusion potentials placed on it by 
the turbulent air. 

In the discussion of the relative merits of the continuum vs. the 
Lagrangian-Eulerian frameworks for developing theories of multiphase flows, one 
of the shortcomings of the continuum approach was immediately pointed out. 
Specifically, the continuum approach requires a priori acceptance of the 
assumption that the particles are of a single size category and have identical 
mechanical properties. Thus, if the particles or droplets are a variety of 
different sizes or vary in mechanical properties then the ’two-fluid’ model has 
become the ’many-fluid’ model. This results in losing one of the original 
benefits of the continuum theories for two-phase flow, that being computational 
efficiency. This raised the question as to whether or not the continuum theories 
are an artifact of the pre-supercomputer era. A conclusion on this question was 
not reached. In defense of the continuum approach to modeling two-phase flow 
it was pointed out that a well posed continuum theory does not contain the 
implicit geometric dependencies that a Lagrangian-Eulerian scheme contains. In 
addition, continuum theories have been used successfully to provide dispersed 
phase field variable information in situations where the added complexity of 
tracking individual particles was not warranted. 
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Due to the stochastic nature and high frequency of collision between particles or 
droplets, a consensus was reached that continuum models may provide a more 
tractable approach than Lagrangian-Eulerian schemes for modeling dense two- 
phase flows. However, advances within this area will require constitutive 
relationships for the dispersed phase spherical and deviatoric stress, separate from 
those of the continuous fluid phase. The development of these relationships will 
have many of the same problems faced in trying to derive viscosity and strain - 
displacement relationships from kinetic gas theory. Advantages of the intuitively 
satisfying Lagrangian-Eulerian approach to multiphase flow theory were discussed. 
The ’particle tracking* allows for the coupling of the fluid and particle at a level 
consistent with single particle hydrodynamic analysis. This approach, despite its 
requirement for increased computational resources can provide a very detailed 
level of information on field variables, including the spectrum of a stochastic 
variable at any point in the domain. Since it is usually the average behavior of 
any given variable, for example particle velocity, that is desired this stochastic 
response is averaged statistically over an ensemble of realizations. The result is 
to produce information at the same level as is available from continuum theory, 
that is, averaged fields. This highlights one of the primary differences between 
the continuum vs. the Lagrangian-Eulerian approach to multiphase flow modeling. 
Continuum theories give averaged field information from averaged conservation 
equations. On the other hand, Lagrangian-Eulerian theories will provide for the 
calculation of stochastic fields, from an ensemble of ’exact’ conservation 
equations. Averaged information from Lagrangian-Eulerian theories requires the 
averaging of fields. 

The modification or modulation of the turbulence structure of the 
continuous fluid phase due to the presence of the dispersed phase was discussed. 
Even at dilute concentrations, the coupling of the phases via drag will result in a 
modification of the production and dissipation of turbulent kinetic energy. 
Questions as to whether or not the presence of the dispersed phase will always 
enhance dissipation of fluid phase turbulent kinetic energy ensued and references 
to experimental observations were made. This led to the discussion of the 
models for the turbulent time scales of the dispersed phase and the continuous 
phase. A passive dispersed phase would have a relaxation time that approaches 
zero and it would follow the continuous phase turbulent and mean motions 
exactly. On the other hand, the relaxation time of the dispersed phase may be 
very large in which case the motion of the dispersed phase would be unaffected 
by the turbulence of the continuous phase. Unfortunately, in many turbulent 
multiphase flows of practical interest there exists a spectrum of dispersed phase 
time and length scales. Further complexity is introduced when it is pointed out 
that, despite the abundant use of single time and length scale models, single 
phase turbulent fluid flow is also dominated by a variety of length and time 
scales over the domain of the flow. Hence, the need for improved models of 
turbulence, for use in theories of single phase flow as well as multiphase flow 
was recognized. Concerning an earlier argument, it was pointed out that regard- 
less of the level of physical detail or lack thereof of the continuum or the 
Lagrangian-Eulerian theory, they both suffer from inadequate turbulence models. 
Any advantage in accuracy gained from using one theoretical approach over the 
other may, in some cases, be washed-out by the inaccuracies inherent in the 
common turbulence model. 

The wisdom and desirability of repeating this workshop at sometime in 
the future was noted in closing. 
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NONGRADIENT DIFFUSION IN PREMIXED TURBULENT FLAMES 


Paul A. Libby 

University of California San Diego 
La Jolla, California 92093 


Abstract - We review recent theoretical and experimental results demonstrating the interac- 
tion between force fields and density inhomogeneities as they arise in premixed turbulent 
flames. In such flames the density fluctuates between two levels, the high density in reactants 
p r and the low density in products p p with the ratio p r /p p on the order of five to ten in flows 
of applied interest. The force fields in such flames arise from the mean pressure drop across 
the flame or from the Reynolds shear stresses in tangential flames with constrained stream- 
lines. The consequence of the interaction is nongradient turbulent transport, countergradient 
in the direction normal to the flame and nongradient in the tangential direction. The theoreti- 
cal basis for these results, the presently available experimental support therefor and the impli- 
cations for other variable density turbulent flows are discussed. 
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INTRODUCTION 


The phenomenology of turbulent flows is largely based on gradient transport assumptions proposed over a 
century ago by de St. Venant and/or Boussinesq. Even in the recently exploited second moment methods 
of analysis, methods which involve large systems of partial differential equations calling for considerable 
computational effort, closure is achieved by employing gradient models to eliminate third moment and 
other quantities. The extension to turbulent flows with variable density of the various gradient models 
carefully developed and validated for constant density turbulence has unfortunately been casually under- 
taken with the consequence that much additional research is needed before the phenomenology is well 
founded for high speed turbulent boundary layers and jets and for turbulent reacting flows, both high and 
low speed. That new significant processes may be operative in turbulent flows with variable density is 
suggested by the recent findings in premixed turbulent flames. It is our purpose to review these findings 
and to discuss their implications. 

Exposition is facilitated if we consider a normal premixed turbulent flame as shown schematically 
in Fig. 1 . Cold reactants, a metastable mixture of fuel and oxidizer, with a mean velocity u 0 , a density 
p r , an intensity of velocity fluctuations characterized by a turbulent kinetic energy k Q and a length scale 
of the large eddies l a are consumed within the flame and exit as hot products with these quantities 
changes to (1 + x) u a , p p , k„ and /„ respectively. Here x is a heat release parameter with values of practi- 
cal interest from five to ten. Figure 1 may be considered an idealization of the flames which occur in 
internal combustion engines, various propulsion systems involving prevaporized fuels and industrial 
accidents. 

The object of a theory for such flames is the prediction of the turbulent flame speed u a and of the 
flame structure. The development of such a theory involves two related considerations; one concerns the 
thermochemistry of the flow, the description of the state of the gas, its density, temperature and composi- 
tion, while the second relates to fluid mechanics. These two aspects are the aerothermochemistry of the 
flow. 
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Various parameters can be used to characterize these flames. One system is proposed by Abraham 
a al. (1985) and involves a Damkohler number Da A = {Mu') (! S L lS L ) where A is a large eddy scale, i.e., 
l o introduced earlier, u' - (3/2 k 0 ) m is a characteristic intensity of the velocity fluctuations and S L and 
5 l are the speed and thickness respectively of an unstrained laminar flame in the chemical system under 
consideration. A second parameter is a turbulent Reynolds number R A = u'Mv where v is a representa- 
tive kinematic viscosity.* A pair of values, Da A and R A , determines the ratio of two velocities u'/S L and 
of two lengths L K I 5, where L K is a characteristic Kolomogoroff length. For premixed turbulent flames 
within internal combustion engines Abraham et al. (1985) determine that 1 < {L K I 8 t ) < 10 2 The impli- 
cation of this finding is that chemical reaction at the molecular level takes place in thin surfaces, laminar 
flamelets, whose structure is dominated by laminar transport and whose motion is determined by the tur- 
bulent velocity field. Similar considerations applied to other premixed turbulent flames indicate that in 
many, indeed most, applications of applied interest this laminar flamelet descnption prevails. 

The Bray-Moss Model and Some Consequences 

Although not restricted to flows involving such flamelets, the Bray-Moss (BM) model simply and with an 
accuracy suitable for many purposes describes the thermochemistry of premixed turbulent flows (cf. Bray 

and Moss 1975). According to this description the instantaneous value of a progress vanable c{x,t ) 

which has the values zero in reactants, unity in products and intermediate values in gas possibly undergo- 
ing chemical reaction determines the entire state of the gas. Thus, for example, 

JL = — 1 — X = 1+tc 0) 

p r 1 + x c T r 

where x is most readily determined from the specialization of the second of these equations to products, 
i.e., from 

T p , (2) 

— = 1 +T 

T r 


* It should be kept in mind that V can vary by one to two 
for determining R A . 


decades within a flame and thus that there is some ambiguity as to its appropriate value 



In averaging the equations for variable density turbulent flows it is useful to employ Favre or mass 
averaging; to illustrate the notion and notation involved consider the i-th velocity component, namely 

pw, 

U,' (x,t ) = = (x) + U (x,t ) = U (x) +M "(x,t ) 

P ~ ~ 

All vanables except the density and pressure are averaged in this fashion. The great advantage of Favre 
averaging is that the conservation equations for variable density flows closely resemble those for constant 
density with exceptions which, as we shall see, indicate special effects associated with the variability of 
the density. Despite this advantage it should not be assumed without support from experiment that the 
closure methods applicable to constant density flows can be carried over to variable density cases without 
modification. Furthermore, although significant only when low turbulence Reynolds numbers must be 
taken into account, Favre averaging is not adaptable without approximation to the molecular terms in the 
conservation equations. The definitive discussion of Favre averaging applied to turbulent combustion is 
given by Jones (1980) while such averaging is emphasized in the monograph on turbulent reacting flows 
by Libby and Williams (1980). Rubesin and cowoikers at the Ames Research Center of NASA apply 
Favre averaging for the analysis of high speed turbulent boundary layers (cf. Wilcox and Rubesin 1980). 

When Favre averaging is applied to Eqs. (1), there results 


JL_ 1 

Pr 1 + X C 



XS 


( 2 ) 


A further important feature of the BM model is an approximation for the probability density func- 
tion (pdf) of the progress variable. In general for a statistically stationary flow 

P (c ;x) = a(x) 5(c) + p(x) 5(1 -c) + tfx) / ( c ;x) (3) 

Now if 

J dc f (c -,x) = 1 

(T 

i.e., if the pdf describing the interior values of the progress variable is normalized, then 
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a(x) + P(x) + ^x)=\ 


(4) 



The essential approximation in the BM-model is that y« 1, an approximation which is satisfied when 
chemical reaction occurs in laminar flamelets as well as in other flow structures, e.g., reaction zones with 
5 l> l k- The implication of this approximation is that the temperature as measured by an idealized ther- 
mometer within a premixed turbulent flame possesses essentially two values, T r within reactants and T p 
within products, and that flamelet passages cause the switch from one level to the other.* 

An important consequence of the y « 1 approximation is exposed if c (x) is calculated from Eq. (3). 
There results 


a(5) = 


1-c 
1 +TC 


1 + X C 


( 5 ) 


where Eq. (4) is used to calculate a(x). We thus see that the strengths of the delta functions is simply 

related to the mean value of the progress variable which must be considered one of the principal depen- 
dent variables in a theory of premixed turbulent flames. 


For the present discussion it is useful to consider the intensity of the density fluctuations; from 
Eqs. (1), (3) and (5) we can compute 



We thus see that the maximum in the relative intensity of the density fluctuations occurs within the flame 
where c = 1/2 and has a value x 2 /4(l + x) = x/4 for x » 1. Thus for degrees of heat release of applied 
interest we must expect intense density fluctuations. 

Although not essential for the present discussion, it is worth noting that the mean rate of creation of 
product w(x) in the BM model involves the product y w max where is the maximum rate of creation 

of product, a quantity determined by the chemical kinetics of the system under consideration. In the BM 
model this product involves a vanishingly small term multiplied by an indefinitely large term and is thus 

* Ex J enSIVe studies the statistics of two valued functions with specific reference to premised turbulent flames have been earned out in older to 
develop a model for the mean rate of creation of product (cf. Bray el al. 1 987 ). 
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indeterminant. As a consequence a separate model for w(x) is needed; progress in this direction is 
described in Bray et al. (1987). 

The Bray-Moss-Libby Model for the Aerothermochemistry 

The notion of bimodality extends to the the description of the velocity components in premixed turbulent 
combustion and leads to the Bray-Moss-Libby (BML) model of the aerothermochemistry of such 
combustion. To appreciate this extension consider the bivariate pdf 

p ( Ui ,c;x) = a(x) 5 (c )P{u l , 0;x) + (3(x) 5(1 - c ) P (u t , l;x) + 

y(x)f(u h c-,x) V) 

Figure 2 shows schematically such a joint pdf. If / («, , c \x) is normalized, Eq. (4) again applies and with 

the approximation y« 1 the statistical behavior of u, and c are again dominated by contributions from 
reactants and products. Note that two of the three pdfs on the right side of Eq. (7) are conditional and 
describe the velocity within reactants and products. It is possible with current diagnostic techniques to 
measure with good spatial resolution the temperature and one or more velocity components and thus to 
determine experimentally the conditional pdfs appearing in Eq. (7) (cf. Moss 1980 and Shepherd and 
Cheng 1988). 

Equation (7) leads to significant results if y 1. To illustrate consider the mean unconditional 
velocity component m, (x) which is given by Eqs. (5) and (7) as 

a i (x) = (\-d)u ir +5 u tp (8) 

where u a and u ip are the conditional mean values of the i-th velocity component within reactants and 
products respectively. We show these values in Fig. 2. Precise definitions of these quantities are as fol- 
lows: 

oo 

u ir (x) = J du, Ui P (m, , 0;x) 
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( 9 ) 

Cj) — J du i P (u ^ , 1 \x) 

— oo 

According to the BM model the mean turbulent flux of all state variables in the i-th coordinate 
direction can be determined from the corresponding flux of the progress variable but from Eq. (7) we 
have 

P u"c" 

— = — = e 0 - c ) ( Uip - u ir ) (10) 

Now a gradient model for this mean flux yields 

P«."c" dc 

p — Vr sr <n) 

where v T is a positive turbulent exchange coefficient. If in connection with Fig. 1 we let correspond to 
x, then dc/dx > 0 implying that p u ”c " < 0 and if gradient transport applies, then u p < u r . But we know 
from our earlier discussion that the velocity in the product stream downstream of the flame where c = 1 is 
(1 + x) times greater than that in the reactant stream. Given this overall increase in mean velocities it 
does not seem reasonable to expect the conditional velocities in products within the flame to be less than 
those within reactants. Accordingly, we see a potential difficulty with a gradient model for premixed tur- 
bulent flames if, as is the case in flames of practical interest, x » 1. 

In Bray et al. (1981), Masuya and Libby (1981) and Libby (1985) the BML model with all gradient 
assumptions avoided is applied to normal and oblique turbulent flames in premixed systems in order to 
clarify the applicability of the gradient model therein. In the next sections we sketch the analyses 
involved and review the results from them 
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APPLICATION TO NORMAL PREMIXED TURBULENT FLAMES 

We now consider application of the notions advanced in the previous section to the flame of Fig. 1. The 
equations needed for the analysis are as follows: 


-j- (pfl) = 0 

dx 




d_,z„ 
dx 


(p u c + p« "c ) = w 


( 12 ) 


and 


d_ 

dx 

d_ 

dx 


pu 

pu 




+ 2 p « " 2 = — 2 u" — Xu 

dx 


pu "c " 


=^- + p u" 2 c" 


du 


dx 
«2 d£ 


+ p«"c"— +p« ^ 


+u"w - 


dx 


(13) 


The usual notation applies to Eqs. (12) and (13) which represent the first and second moment equations 
respectively for one dimensional flames. The quantities y> u and describe dissipation effects which for 
the flames under consideration are due solely to chemical reaction, i.e., are proportional to w . Similarly 
the velocity-chemical reaction term u 77 ^ can be convincingly described by application of the BML 
model; it is found to be proportional to w, models for which are presently under development (cf. Bray et 
al. 1987) but which are not needed for present purposes.* 


Some analysis permits Eqs. (12) and (13) to be reduced to two equations with c as the independent 
variable and with a dimensionless velocity intensity I =pu" 2 /p r u 2 and a dimensionless turbulent flux 
F = pu "c "lp r U Q as the two dependent variables. In the present discussion the latter variable is of more 
significance. In this formulation w is needed only if subsequent to finding the solutions for 1(c) and 
F(c) the spatial distributions are sought. However, comparison between theory and experiment is 

. It is worth noting that in these equations the effects of pressure fluctuations are neglected. In constant density turbulence models for such ef- 
fects have been painfully and carefully developed over a period of many years but these models do not apply to reacting flows with heat re ease 
!mte ^bt.rmduced pressure fluctuations result from volumetric sources, i.e.. from an entirely different mechamsm. Appropnate models 

for such fluctuations are not presently available. 
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conveniently carried out in terms of c so that generally the absence of a model for w is not a shortcoming 
relative to the analysis of one dimensional flames. For premixed turbulent flows in two and three dimen- 
sions such a model is needed and research to that end is underway (cf. Bray et al. 1987). 


In developing the final equations for /(c) and F (c) the first of Eqs. (2) and (12) permit u to be 
eliminated while the second and third of Eqs. (12) are used to eliminate dp/dx and vv respectively. Clo- 
sure requires the third moment quantities p m" 3 and p u" 2 c", u", c" and the dissipation terms to be 
expressed in terms of the two dependent and the independent variables but with the exception of some 
inessential uncertainties in the conditional velocity statistics suitable models for these quantities free of 
gradient assumptions can be developed. For example, we have 




C" = T C(1 — 
1 + XC 


(14) 


In the present discussion two terms in Eqs. (13) are of particular interest; we refer to those involving 
dpldx. If Eqs. (12) and (13) are specialized to constant density flows, e.g., by setting x = 0, then from 
Eqs. (14) we see that these terms vanish and indeed the reduction to standard equations for such flows is 
complete. The implication to be drawn from this consideration is that these pressure gradient terms 
account for effects operative in variable density turbulence, in particular an interaction between a force 
field associated with the pressure gradient and density fluctuations associated with heat release. In 
premixed turbulent flames this interaction is present only within the flame proper since both the force 
field and the density fluctuations are absent in the reactant and product streams on each side of the flame. 

Additional comments on this interaction are indicated. If conventional rather than Favre averaging 
is used, the interaction is contained within the resulting equations but is obscured by the clutter of a large 
number of terms involving density fluctuations. In more general variable density flows models for the 
multipliers of the several pressure gradient terms must be introduced as indicated for nonpremixed 
combustion by Jones (1980) and for high speed flows by Wilcox and Rubesin (1980). 


The final equations for 1(c) and F(c) involve singularities at c =0,1 as is to be expected since 
these points correspond to points at infinity in the x-variable but appropriate series expansions permit the 
numerical solutions to be initiated in the neighborhood of these end points. A shortcoming of the theory 
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is that I a = u' 2 a /u 2 where u' 2 a is the intensity in the reactant stream upstream of the flame must be 
imposed, i.e., the turbulent flame speed is not calculated. This shortcoming can be turned into a virtue for 
the purpose of studying the structure of premixed turbulent flames since it can be argued that selection of 
I a so as to achieve agreement with experimental results assures that the predicted flame structure 
corresponds to realistic flames. In our original work we set I a = 0.22 unless we had specific reasons to do 
otherwise but in recent years our confidence in this value has been reduced by evidence that a wide 
variety of values is found depending on the geometry of the flames and on the flameholding mechanism 
employed. The detailed reasons for this ambiguity are unknown.* 

Numerical Results for Normal Flames 

In Fig. 3 we show the distribution of the turbulent flux in terms of F(S) from both theory and from the 
experiment of Moss (1980). In the open flame studied by Moss x = 6.5 and I a =0.16. The two curves 
represent theoretical predictions based on slightly different models for the conditional statistics, differ- 
ences which are irrelevant for the present discussion. If it is recalled that gradient transport for these 
flames implies F < 0, we see that the entire flame structure exhibits countergradient diffusion. The expla- 
nation for such diffusion resides in the interaction between the pressure drop across the flame and the den- 
sity fluctuations which from Eq. (6) are seen to involve a relative intensity greater than unity, roughly 1.4. 
The pressure drop tends to accelerate the light parcels of product relative to the heavy parcels of reactants 
in the direction of high product concentration, a differential effect contrary to gradient transport. 

Several further comments are indicated. Calculations with the pressure gradient terms in Eqs. (13) 
suppressed yield F < 0 for all values of x. Moreover, as x — » 0 the theory predicts gradient transport and 
thus the expected behavior for nearly constant density turbulent flows. These results lend credibility to 
the explanation of interaction as the cause of countergradient diffusion and to the validity of the theory in 
general. The notion of countergradient transport in premixed turbulent flames is now widely accepted 


* A recent study (Libby 1987) removes the shortcoming with respect to I 0 by invoking the Hakberg-Gosman condition and discusses in some de- 
tail the uncertainties in its value. 


146 


and has been observed in a variety of laboratory flames with simple geometries. Later we discuss its 
more general applicability. 
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THE OBLIQUE PREMIXED TURBULENT FLAME WITH CONSTRAINED MEAN STREAM- 
LINES 


In Fig. 4 we show schematically a premixed turbulent flame which is oblique to the reactant stream. Such 
flames can be established by a flameholder, e.g., a cylinder, in a duct carrying reactants. The classic 
experiment of this nature is due to Wright and Zukoski (1962) while the theoretical description of these 
flames based on the BML model is given by Masuya and Libby (1981). If we assume that the duct prohi- 
bits significant deflection of the mean streamlines both upstream and downstream of the flame, it is rea- 
sonable to consider as an idealization an infinite planar flame held at a specified angle 6 with undeflected 
mean streamlines. In this case description of the flow involves an analysis identical with that for purely 
normal flames and a second, subsequent analysis of the tangential velocity component which involves 
explicitly the flame angle 0. The treatment of the tangential flow requires determination of the intensity 
of the fluctuations of the v-velocity component and the tangential flux of the progress variable which in 
dimensionless form are I v /, p r u 0 2 and F v S ^V7p r respectively. This latter variable is of 

principal concern in the context of the present discussion and from Eq.) is given by 


pv c 

=c o -o Op -v r ) 

where v r and v p are the conditional tangential velocity components within 
tively. 


(15) 


reactants and products respec- 


Gradient transport indicates that this mean flux is zero since the tangential gradient of the progress 
variable is zero and thus that the two conditional velocities in the tangential direction are equal. The 
implication from this result is that parcels of reactants and products have streamlines which differ only by 
the differences in the noimal conditional velocities. However, there is a tangential force field in these 
flames arising from the x-wise gradient of the Reynolds shear stress ^?V'; the existence of this forces 
field can be seen from the mean conservation equation for tangential momentum which yields 


dx tan 0 dx 


Note that when the heat release and thus the density fluctuations vanish, this force field is absent. 


( 16 ) 


This discussion establishes that oblique turbulent flames in premixed systems involve a tangential 
force field and density fluctuations. Thus according to our previous argument we can expect nongradient 
turbulent transport in the tangential direction. Indeed calculations show that v p > v r , that pv"c " > 0 and 
that the parcels of reactants and products follow different mean streamlines with the former exhibiting 
only small tangential velocity while the latter possesses large tangential velocities. This behavior is 
shown schematically in Fig. 4. We thus have an example of nongradient transport to add to the previously 
discussed case of countergradient transport. To date there have been no detailed measurements in oblique 
flames to assess the validity of this theory. 
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CONCLUDING REMARKS 


We show that when mean force fields and density fluctuations coexist, an interaction leads to turbulent 
transport which is not described by the usual gradient model. In the simple flow configurations associated 
with normal and oblique planar flames, the latter with constrained mean streamlines, countergradient and 
nongradient transport exists. There is no question that the notions suggested by these findings are con- 
ceptually important. Moreover, the experimental results of Heitor et al. (1987) relative to the complex 
flow associated with a baffle stabilized premixed turbulent flame establish that nongradient transport 
exists in complex flow configurations, i.e., in flows of applied interest. Within the context of the present 
discussion the first few sentences of this paper are worth quoting: 

In turbulent, premixed flames there arise source terms, explicitly set out below, in the conservation 
equations for the turbulent heat transfer rate and stresses that have no counterparts in non-reacting 
flows. Analysis (Masuya and Libby 1981; Bray, Libby and Moss 1985; Libby 1985) has shown that 
at least in the two idealized extremes with the flame either normal or oblique to the approaching 
reactants, and at practically important levels of heat release, these terms are sufficiently large to 
cause non-gradient transport of turbulent heat flux. This finding is important because it casts doubt 
on the applicability of turbulence models that use gradient-transport hypothesis. 

From a detailed study of the velocity and temperature fields Heitor et al. conclude that the interaction 
between pressure gradients and density fluctuations results in the largest contribution to the balance of the 
turbulent heat flux and that that flux is not aligned with the mean temperature gradient 

The importance of nongradient transport on general variable density turbulence remains to be esta- 
blished. With respect to premixed turbulent flames we note that the question arises as to the influence of 
these processes arising within their structure on their orientation between reactant and product streams 
remains to be clarified. Answering this question requires further theoretical and computational efforts.* 

* If the provisional model for W (x) set forth in Bray el al. (1987) is validated, the effort required is largely computational. 
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Additional research is also required to determine the importance of the interaction under discussion 
in high speed flows. It would seem obvious that in scramjets with their large gradients of mean pressure 
and mean turbulent shear stresses this interaction must be operative but whether it plays an important role 
in the mixing and chemical processes in the flow is unclear. Certainly the present casual, uncritical exten- 
sion to supersonic chemically reacting flows of the phenomenology of constant density turbulence should 
raise healthy skepticism concerning the validity of the resulting predictions but unfortunately that skepti- 
cism does not appear to be shared by authors, reviewers and editors of current journal articles.** Similar 
uncertainties would appear to prevail for turbulent boundary layers in high speed flows so that additional 
research is indicated. 


** Occasionally there is a mild rejoinder warning the reader of potential fundamental difficulties. An example from a recent report is as follows: 
The turbulence model we we have used is the fc — £ model as described in ... Although Reference ... specifically addresses the question of 
compressible flows, it should be noted that turbulence models for compressible flows are not well-developed. Hence, the choice of the standard 
k — £ model in this situation cannot be regarded as definite." 
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THE TURBULENT HEAT FLUX IN LOW MACH NUMBER FLOWS 
WITH LARGE DENSITY VARIATIONS 


Peter J. O’Rourke and Lance R. Collins 
Theoretical Division, Group T-3 
Los Alamos National Laboratory 
Los Alamos, New Mexico 87545 


L INTRODUCTION: THE DIRECTED ENERGY FLUX 

fnr f J h l S JHP er iS r? n f r . ned W . ith a Physical effect of fundamental importance 
Thp in ^ tur ^ u ^ nc ? transport in flows with large density variations. 

•^ e f ffect occurs because the interaction of pressure and density gradients gives 
rise to a turbulent heat flux, which we call a directed flux, that is not accounted 
for in turbulence models for constant density flows. To see how this flux arises it 
is perhaps best to consider an example of Rayleigh-Taylor instability, as depicted 
" ^ , g ’ A h ea yy, cold gas overlays a light, hot gas in a box with a gravitational 

pnt !/ | 10n / n a< r gatlve z ". direc ti°n. The induced hydrostatic pressure gradi 
mitin^ e T^ ateS i th ® J‘ ght gas in ]° the heav y gas and causes the instability and 
» n A T “ e ™ ]oc ' tles averaged across a horizontal plane are in the z-direction 
mi vi h® cau a e the heavy gases are falling, the mass-averaged velocities in the 
ing region will be negative. Relative to a surface moving downward with the 
mass-averaged velocity, there will be a net upward flux of energy. This is 

thp^L a n th °Tk the maSS ° f L ight gas crossin g th e surface upward equals 
the mass flux ofheavy gas crossing downward, the light gas, being hotter, carries 
with it more energy per unit mass. B s 

This upward energy flux is the directed flux. In the example of Fig. 1 this 
Ji m d i r rectl . on o f the negative of the temperature gradient, just as 
Srn th? la mmar Fourier heat conduction law. If after some time we were to 

thp twn o-oc T r /° ? a . 4 th J lght g ? s ove , rla y the heavy gas, to the extent that 
t e two gases had not already mixed on the molecular level, there would be an 

unmixing in which the light gas would separate from the heavy gas. In this 





Z 


gravity 

Fig. 1. Schematic depiction of 
Rayleigh-Taylor instability. 
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unmixing process, the directed heat flux would be upward, in the direction of the 
mean temperature gradient and opposite the direction given by the Fourier heat 
conduction law. This phenomenon has thus been called countergradient delu- 
sion. 1 It cannot be predicted by turbulence models that use gradient transport, or 
a Fourier-like law, to describe turbulent heat transport. In our example, the tur- 
bulent heat flux <t>^ is in the direction opposed to the pressure gradient, rather 
than the temperature gradient. We shall see that taking «(>« ~ — Vp is often more 

realistic for turbulent flows. . , „ 

Between single-phase, two-density turbulent flows and two-phase Hows 
there is an analogy that we will exploit in our turbulence modeling. This analogy 
will be used in helping to formulate the equations and in developing closure 
approximations for some of the terms. In two-phase flow modeling, which has 
received much attention within the last ten years2,3 separate mass, momentum 
and energy equations are kept for each phase, and these equations are coupled 
through functions that give the exchange rates of mass, momentum and energy 
between the phases. Following Besnard, Harlow, and Rauenzahn,4 we will use 
an alternative, analogous formulation. In place of two mass equations, we will 
keep an equation for the mean density and one for density fluctuations. In place 
of two momentum equations, we will keep a mean momentum equation and an 
equation for mean velocity differences associated with fluid elements of dinering 
density. We use this second formulation because it allows for the possibility ot 
modeling not just two-density flows, but the flows with a spectrum of densities 
that often occur in practical applications. . n 

Two physical examples, one of two-phase flow and one of single-phase flow 
with density variations, serve to illustrate the analogy and another situation in 
which countergradient transport can arise. In both examples, pressure gradients 
are responsible for centrifuging lighter material inward toward the centers ot 
rotating flows. In a two-phase bubbly flow, this effect has been observed in the 
vortices in the wake of an obstacle. 5 In single-phase flow, it is probably responsi- 
ble for the experimental results of Wahiduzzaman and Ferguson.6 The experi- 
menters measured the radial temperature profiles in an axisymmetric swirling 
flow in a constant volume cylinder. The experimentally measured temperatures 
are plotted at four different times as the circles in Fig. 2. The lines are computed 
temperature profiles using the KIVA code? with a fe-e turbulence model» and 
gradient heat transport with a turbulent Prandtl number of 0.9. It can be seen 
that a hot region in the center of the cylinder persists much longer in the experi- 
ment than in the calculation, showing the large errors that can arise when a 
gradient heat transport approximation is used. . 

The phenomenon of countergradient transport in single-phase flows was 
recognized seven years ago in research on the structure of turbulent premixed 
flames l In retrospect, it is easy to see how this phenomenon arises. Figure 3 
depicts schematically a planar turbulent premixed flame with velocities shown in 
the frame of reference of the flame. Mass conservation and the fact that the corm 
bustion is nearly isobaric together imply that the hot product gas velocities will 
be larger than those in the reactants. Momentum conservation then implies that 
the pressure in the products will be lower than in the reactants. Since the dir- 
ected heat flux is in the direction opposed to the pressure gradient, this heat flux 
will be from the colder to the hotter gases; that is, it will be countergradxent 

Two approaches have been used for modeling turbulent premixed flames -- a 
single-phase formulation and a two-phase formulation. In the single-phase for- 
mulation of Bray, Moss, and Libby (BML formulation), equations are kept 
for the mean product gas concentration, the mean momentum, the turbulent 
fluxes of these quantities, and for the dissipation rate of turbulent kinetic energy. 
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Fig. 3. Schematic depiction of a turbulent premixed flame. 

Spalding 13 utilizes a two-phase formulation, retaining mass, momentum, and 
energy equations for each phase. Spalding 13 assumes that the major source of 
mixing is due to the difference between the mean velocity of the phases and thus 
ignores the turbulent kinetic energy within each phase. The BML formulation 
accounts for both sources to the turbulent kinetic energy. In principle, the equa- 
tions of one formulation should be derivable in terms of those of the other, al- 
though to our knowledge such a derivation and comparison have not been made. 
In this paper some closure approximations are proposed, based on a derivation of 
the single-phase equations from two-phase equations. Only the BML formulation 
has been compared with experimental measurements of turbulent flames, 1 1 and 
satisfactory agreement was obtained. 

In practical applications of turbulent combustion, other physical effects that 
cause mixing and unmixing are superimposed on the pressure drop across the 
flame. At Los Alamos, we have been involved for the past twelve years in the 
numerical modeling of combustion in internal combustion engines. 7 * 13 - 16 
Figure 4 illustrates some of the complexities of the turbulence/chemistry inter- 
action in an engine burning premixed charges. A turbulent premixed name is 
propagating away from an ignitor located near the center of the cylinder head 
wall. Mach numbers are small, and thus the mean pressure is nearly uniform in 
space 17 and changing with time due to piston motion, combustion, and wall heat 
loss. Near the top of its motion, the piston, and the axial flow velocities in the 
combustion chamber, decelerate. This causes a small positive axial pressure 
gradient and induces Rayleigh-Taylor instability and mixing where the flame is 
propagating downward in the axial direction. This same pressure gradient will 
cause a differential axial acceleration of the hot products and cool reactants and 
promote Kelvin-Helmholtz instability where the flame is propagating radially. 
Swirl, a nearly symmetric rotational motion of the burning gases, is introduced 
by engine designers to promote mixing but will have two competing effects in the 
engine of Fig. 4. Swirl induced shears will enhance turbulence and mixing, but 
the radial pressure gradient caused by the swirl will, as in the experiments of 
Wahiduzzaman and Ferguson, 6 cause countergradient transport and suppress 
mixing. It is important to point out that among these various turbulence effects, 
only those associated with shear instability and mixing are accounted for in 
current engine models. 

In predicting turbulence in internal combustion engines and other practical 
combustors, one cannot use two-phase models or single-phase equations for two 
density flows. This is because within the reactants and products there will be 
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cylinder head 



Fig. 4. The turbulence/chemistry interaction in an internal combustion engine. 

distributions of density. In the engine of Fig. 4 density differences will arise due 
to wall heat transfer and due to entropy differences in the product gases of this 
confined burn. 18 In other combustors density distributions are caused by charge 
non-uniformity and spray vaporization. 

In the next section we derive preliminary equations for a turbulent fluid 
with large density variations. Our aim is to develop a model that has three 
attributes: 

(1) the model can predict mixing and unmixing due to shear instabilities and 
pressure-density gradient interactions; 

(2) the model can account for a distribution of densities; and 

(3) the model equations can be efficiently integrated in two and three 
dimensions. 

The second attribute precludes use of two-phase flow equations, although investi- 
gation of the two-density limit will yield valuable information. The third attri- 
bute precludes use of the Reynolds stress equations, especially in three dimen- 
sions. It seems appropriate to seek a simple one- or two-equation extension of 
popular two-equation models for turbulent shear flows. 

II. THE EQUATIONS 

A. Overview 

We first derive equations for the average density p and the Favre-averaged 
velocity u and enthalpy h. Using the low Mach number assumption, we relate 
the turbulent heat flux to the difference between the average velocity u and 
the Favre-averaged velocity u. We denote this difference u — u by a, and the 
transport equation for a is derived and discussed. Closure approximations for 
terms in the a-equation are postulated based on the analogy between two-phase 
flows and single-phase, two-density flows. A comparison between the single- 
phase and two-phase equations suggests that the fluctuating stress terms in the 
a-equation are primarily associated with the decay of a. We present an algebraic 
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closure approximation for a that results in a heat flux that is the sum of contribu 
tions due to gradient and directed transport. Our a-equation is compared with 
similar equations in the literature. 

B. The Equations of a Low Mach Number Flow with Large Den sity Variations 
For the low Mach number flow of a single component ideal gas with large 
density variations, the equations are the following: 


ap 

— + V • (pu) = 0 , 

dt 

(continuity) 

(i) 

+ V-(puu) + Vp = V- o + pg , 
dt 

(momentum) 

(2) 

ap h , dP 

h V • (pu h) = — — 1- Q , 

dt dt 

(enthalpy) 

(3) 

P (t) M ~ pRT , 

IV 

(thermal equation of state) 

(4) 

h(T) = [ cWdi, 

(caloric equation of state) 

(5) 


where o is the laminar viscous stress tensor, Q is the volumetric rate of heating 
due to such sources as chemical reaction or divergence of the laminar heat flux, 
P(t) is the volume average pressure of the system, p(x,f) is the pressure fluctua- 
tion from the mean value P(t), and M w is the molecular weight of the gas. For low 
Mach number flows Ipl/P = M 2 , where M is the Mach number, 17 From these 
equations one can derive an equation for the divergence of the velocity field: 


Vu 


l dP 

yP dt 



( 6 ) 


In an open flow system, P is just the ambient pressure; for flow in a closed volume 
V, an equation for P can be derived by integrating (6) over V: 


= + I 1- Qdu . 

yP dt V dt yP V J v 

C. The Averaged Equations . , , , ,. 

In our turbulence equations we will use both unweighted averaged quanti- 
ties and Favre (density-weighted) averaged quantities. The unweighted average 
and Favre average of a quantity 4> are defined respectively by 


- l V' 

<b (x , t) = — > <p (x,f) 
^ NE “ 


( 8 ) 


and 
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1 


(9) 


0(x , t) - 


P (x , t) 


1 

NK 


X P a (x,t)<p a (x,t) 


where ensemble averages are used, <p a being the value of 0 in a particular 
experiment a and NE being the number of experiments in the ensemble. The 
fluctuations from these averages are denoted by 


4>' = <P a - <P 


( 10 ) 


and 


<t>" = <t> a -4> , ( 11 ) 

where we drop the subscript a on the fluctuations. 

By averaging Eqs. (l)-(5) we obtain the turbulence equations: 


7+MpuVo, (12) 

OL 


dp 11 

dt 


+ V • ( p uu) + Vp = T • K + V • a + pg. 


(13) 


dp h d P - . 

— + ▼•( P M= — + Q - V-4> h , (14) 


PM = pRT , 

IV 


(15) 


and 


h = I c{ i) dx , (16) 

where th e Reyn olds stress R is given by - pu'u' and the turbulent heat flux <1>* is 
given by pITTF! In deriving Eq. (16) we have assumed that the characteristic tem- 
peratures over which significant changes in c p occur, are much larger than char- 
acteristic temperature fluctuations T \ 

D. An Alternative Expression for <b^ 

~ An alternative expression for the turbulent heat flux can be derived from 
the averaged and unaveraged equations of state. Subtracting (15) from (4) 
results in 


P' M w — pRT — p RT = pRT" + p'RT . (17) 

We now assume that 


|P'|/ P << ip'l/ p . (18) 
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Under this assumption the left-hand side of (17) can be neglected, and we obtain 



(19) 


wherein temperature and density fluctuations are directly related. 

The assumption (18) requires some discussion. It will certainly be true for 
an open flow system, because the P is always the ambient pressure and never 
fluctuates. In a closed system, such as in internal combustion engine cylinder, 
there can be significant fluctuations in the mean pressure P, primarily due to 
fluctuations in the chemical heat release rate. In an internal combustion engine 
these are referred to as cycle-to-cycle variations, and there is currently some 
debate 20 whether or not these cycle-to-cycle variations should be called turbu- 
lence. We recommend that in performing the averages (8) and (9) one should use 
only experiments for which the mean pressure history P(t) is nearly equal to P(t). 
When this limited ensemble average is used, the assumption (18) is automat 
ically satisfied. 

Subtracting (16) from (5) gives 


f T+T ‘ ~ /on' 

h'= L c(x)di=c(T)T\ 

Jr p p 

where again we use the assumption that c p varies little for temperature changes 
equal to T. Using (20) and (19) in the defining formula for gives 


= pu "h" 


= -c(T)T p’u" 

= pclflTa 
p 


( 21 ) 


where 


- - ~ ( 22 ) 

a = u = u — u . 

Equation (21) is the alternative expression we seek for the heat flux. It says 
that <|>Si is proportional to a quantity a that can be loosely thought of as the differ- 
ence between the volume-averaged and mass-averaged velocities. In a two- 
density or two-phase system, a is proportional to the difference between the veloc- 
ities of the two phases. In order to investigate further the nature of the turbulent 
heat flux, we must derive a transport equation for a. 

E. The Transport Equation for a . 

The transport equation for a is obtained by subtracting the equation tor the 
mass-averaged velocity u from that for the volume-averaged velocity u. The 
result is 
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da ~ ~ / 1 1 \ - 1 

f u-Vu — u • Vu + i - — 3 ) V p + — Vp' 

d t \ P p / p 


= - =V*R +(- - iW a + - V-o' . (23) 

p V P p / P 

As an aid in modeling some of the terms in Eq. (23) we will derive an 
a-equation for two-phase flows and compare this with Eq. (23). It will turn out 
that all of the terms in Eq. (23), with the exception of those associated with the 
fluctuating stresses p' and o', will be d uplicat ed e xactly b y terms in the 
a-equation for two-phase flows. The terms 1/p Vp' and 1/p V • o' are then associ- 
ated with terms in the two-phase a-equation that arise due to momentum 
exchange between the phases. With this comparison as a guide we postulate a 
model for the fluctuating stress terms. 

We now briefly introduce the equations for two-phase flow. For more details 
the interested reader should consult Refs. 2 and 3. For simplicity we restrict our- 
selves to the flow of two incompressible phases. The continuity equations for each 
phase are 


dp i Q i 

— + V-( Pl a lU ) = J 2X 


(24) 


and 


dp 2 Q 2 

dt 


+ V 


( P 2 a 2 u2) = J 12 


-J 


21 


(25) 


In these equations pi are the microscopic (or conditional) densities, which we 
assume to be constant; ui are the average velocities within each phase (where the 
"t” next to the overbar indicates a conditional average in phase i); a, are the vol- 
ume fractions of each phase (ai + a 2 — 1); and J 21 is the rate of mass transfer per 
unit volume from phase 2 to phase 1. The momentum equations for each phase 
are 


a Pi Q i u -1-1 

h V • (PjQjU u l + a ( Tp = p^jg + OjV • o + V • (a^j) + P 


dt 


(26) 


and 


ap 2 Q 2 u2 - 2-2 

— — + V • (p 2 a 2 u u ) + a 2 V p = p 2 a^g + a 2 V • o + V • (a^ - P 2J . 


(27) 


As is commonly done in two-phase flow modeling, 2 we assume here that the 
phases are in local pressure equilibrium; that is p 1 = p 2 = p. The conditional 
Reynolds stresses Rj are given by 
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— I — l 

R ( = — p ( ( u — u ) (u — u ) . 

The rate of momentum exchange per unit volume from phase 2 to phase 1 is 
denoted by I* 2 i- 

Average flow variables f are related to the averages within each phase f l by 

f = a / + a ,/ 2 

Thus, for example, 

P = Vl +Q 2 P2 


(28) 

(29) 


and 


- ~ —1 —2 
pu — QjPjU + a 2 p 2 u 

By using the relations (29) and adding Eqs. (24) and (25) and Eqs. (26) and (27) we 
obtain the same total mass and momentum equations, Eqs. (12) and (13), that we 
have previously derived, when it is realized that 


K = a 1 K 1 + a 2 K 2 


- P jQ |(u * - u)(u’ - u) 

- p 2 a 2 (^ - u)(u 2 - u) . (30) 

Thus the mass-averaged velocity equation is 


du ~ ~ 1 - 1 - 1 

— + u Vu + = V p = = V • a + = V - R + g . 

dt p p p 


(31) 


To obtain the a-equation, we will subtract (31) from an equation for u, 
which will now be derived. By dividing (26) by pi and (27) by p2 and summing the 
results, we obtain 


du —1—1 —2—2 / 1 \ — 

b V • (a . u u + a.,u u ) + I — ) V p 

dt 1 2 ' p/ 


= g + 



(° 1 a 2 \ 
+ V- ( Rh Rj 

^P, 1 P., l > 


+ p -(7- 



(32) 


It can be seen that 
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~V (uu)= u*Vu + uV u . 


(33) 


- 1-1 - 2-2 G 1 a 2 
a t u u + a u u - — K - — K 


1 


P.> 


Using (33) in (32) and subtracting (31) yields the two-phase a-equation: 


da - ~ / 1 1 \ - 

— + u-Vu -u-Vu + uV*u + - - zV p 

dt \ p p ; 


1 / 1 

= — = V • R + ( - 

p v P 



o + V 


21 




Comparing (34) and (23), we see they agree if 


(34) 


l l / l l \ 

- v P ' - - V o' = uV-u + P... — - — . (35) 

P P V P 2 

This is the relationship we seek between the fluctuating stress terms and the two- 
phase momentum transfer terms. To obtain closure we need to postulate a form 
for 1*21 > and we will investigate expressions employed in two-phase flow 
modeling. 

F. Expressions for the Momentum Exchange Rate in Two-Phase Flows 

For a dispersed phase 2 of equal-sized spherical particles of radius r in a 
continuous phase 1, an expression for the momentum transfer term is 2 1 


I* = - c 

21 8 C l) 


P l °2 2 — 1 


2 — 1 . —2 —1 —2 
U — U I (u — U ) + J 9) U 


(36) 


Thus, P 21 has two terms — one due to aerodynamic drag and one due to mass 
exchange. This form of the mass exchange term assumes there are no circulation 
velocities within the particles. Equation (36) has theoretical justifications 1 when 
P2 > > pi and velocities within each phase are sharply peaked near their mean 
values. It neglects virtual mass effects, Basset history effects, and particle dis- 
tortions and oscillations. 22 

Motivated by Eq. (36), modelers usually use a similar expression for all two- 
phase regimes: 2 




(37) 


where K is called the drag function and u« is some average interface velocity. The 
quantity K is a positive function of pi, p2, ai, a 2 , |u 2 — ufi, and an entity size r. 

If we accept Eq. (37) then one is led to the postulate that the fluctuating 
stress terms in the a-equation (23) contribute to the decay of a. Indeed one can 
show that for a two-phase flow 


a = 


V2 ( P, ~ P 2> 


-2 

(u 


p 


— 1 

u ) , 


(38) 
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and (38), in conjunction with (35) and (37), gives 


l 1 ( 1 l 

- V p' — - V • o' - uV u + K ' a + — - — 

P P 21 V P 2 P x 

where 


(39) 


p 

K = - K . 

V2P1P2 

K has dimensions of a frequency. 

G. Final Form of the a-Equation 

After substituting (39) in (23) and some rearrangement of terms one obtains 

~ ~ / 1 1 \s~~b — 

— + V • (ua + au) + K'a + J — - — u + uV • u + = V p 

dt 1 V P 2 Pi ’ P 

1 1 - 1 

- = V • R - V • (u'iO = = u"u" • V p + = V • ( p'u"u’ ) . (40) 

P p p 

Here we have introduced the quantity b as a dimensionless measure of the 
density fluctuations: 



If the density fluctuations are not too large, then b is approximately a self- 
correlation coefficient for density fluctuations: 


b = 



In fact. Ref. 4 uses 


(42) 


B = (p ') 2 ( 43 ) 

as a measure of the density fluctuations. We will develop a transport equation for 
b in future work. 

Three further terms in (40) must be modeled. First we deal with the mass 
exchange term. One can show from (24) and (25) that 


V 




(44) 


and hence the fourth and fifth terms on the left-hand side of (40) combine to give 
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(45) 


uV • u + J ( — - — )u" = J ( — - — W - u) - uV • a 
21 V P 2 P, > 21 V P 2 P, ' 

We assume that u s = u. An assessment of the validity of this assumption must 
await a precise physical int erpr etation of the quantity u s . 

Second, for the tensor u"u" one can show that for two-phase flows 



where B is defined by (43) and given in two-phase flows by 

R = - p 2 ) 2 . 

For future reference we also note that 


(46) 


(47a) 


6 = — (47b) 
P1P2 

in two-phase flows. We define the volume-averaged conditional Reynolds stress r 

by 


R ! R 2 (48) 

r = a — + o, — 

' P, P‘2 

As a first approximation, and despite experimental evidence to the contrary in 
turbulent flame experiments,® we assume the conditional Reynolds stresses are 
equal and isotropic. Then 


= (49) 
P, P 2 3 

where k' is related to the specific turbulent kinetic energy k = i(u") 2 by 


A transport equation for k will be developed in future work. 

We also use a two-density distribution to model the triple correlation term 
in Eq. (40). After some algebraic manipulation and use of the assumption (49) 
one obtains 


-2 

= + jf)“- 


(51) 
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By substituting (45), (46), (49), and (51) in (40) and using the approximation (42) 
we obtain the final form of the a-equation: 


da ~ ~ ~ u — 

— + u • Va + aV • u + a • Vu + V • (aa) + - V p 
dt p 


aa V p 2 

= -K'a+ — + ~{k 

On 3 


1 — \ X_p 

2 b I p 


(52) 


H. An Algebraic Closure Approximation . 

In numerical computations of multidimensional fluid flows, use ot Eq. (oZ) 
would require solving two or three additional transport equations for components 
of a Although this is not an unrealistic task for modern computers, considerable 
computational efficiency would result if an accurate algebraic closure approxima- 
tion for a could be found. In this section we present such an approximation based 
on an assumption whose validity must be tested in experimental comparisons. 
The resulting expression for a predicts gradient heat transport, but also contains 
a contribution that predicts the directed flux arising from the interaction of pres 
sure gradients and density inhomogeneities. f 

The assumption we make is analogous to the drift flux approximation ot 
two-phase flow modeling. 2 In two-phase modeling this assumption is that the two 
velocity fields are so tightly coupled through the drag terms that characteristic 
drag times are much smaller than characteristic flow times. For us the assump- 
tion is that 


where u 0 and L are a characteristic velocity and gradient length for the flow. 

Assuming (53) is true, order of magnitude estimates of the terms in (52) 
show all terms on the left-hand side can be neglected, except the pressure grad 
ient term. On the right-hand side, the dyadic product term aa is negligible since 
Eqs. (38) and (47) show that a/6 is proportional to the velocity difference between 
fluid elements of different density. The resulting equation for a becomes 


a 


1 

~K' 


b - 2 

= V p + - k 

P 6 



Equation (54) can be put in a more recognizable form if we use 


(54) 


Vp 

p 



a = 


1 

~K' 



(55) 


In conjunction with (21), Eq. (55) gives a heat flux that is the sum of contributions 
proportional to -Vpand -VT. The former is the directed flux. It goes to zero in 
the absence of density fluctuations b. 
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The gradient transport term in (55) looks similar to the gradient heat flux 
commonly used in turbulence modeling, but there is a difference. The usual form 
used for the turbulent heat flux8 is 


= - pc (56) 

p Pr T c 

where Pr,y is the turbulent Prandtl number and t the turbulence dissipation rate. 
Equation (56) agrees with the heat flux contribution obtained from the second 
term in (55) if the drag time associated with fluid elements of differing density 
equals the turbulence dissipation time. For the momentum exchange function 
(36) it can be seen that the drag time is 


1 P2 


r 


K' 


p, 


r» _2 i 

u - u 


(57) 


when c D is approximately unity. On the other hand, the turbulence dissipation 
time is 


k 

E 

where L is a turbulence length scale. Equations (57) and (58) agree if r == L, ~k> = 
lu 1 - u2| and pi ~ p 2 , but when these equalities are violated more accurate heat 
fluxes could be obtained using a drag time, and not a turbulence dissipation time, 
to evaluate the heat flux vector. 

III. COMPARISON WITH OTHER WORK 

In this section we compare our a-equation with two others in the literature. 
In the BML formulation for turbulent flames, 10 an equation is kept for the tur 
bulent flux of reaction progress variable c. Our quantity a is just a constant times 
the turbulent flux of c: 


I. 


(58) 


(59) 


where p r and p p are the reactant and product densities. Two differences are ob- 
served between the a-equation one derives from the BML formulation and ours. 
First, in the BML formulation it is not assumed that the conditional Reynolds 
stresses within each phase are equal and isotropic. An equation for the uncondi- 
tional Reynolds stress R is retained, and the difference between the conditional 
Reynolds stresses is modeled using R. Accordingly, the double and triple correla- 
tion terms on the right-hand side of (40) are modeled in a more detailed fashion 
although the authors observe 10 that "this modeling is generally not found to be 
too critical to the predictions of first- and second-moment unconditional 
quantities.” 

The second difference is in the modeling of the fluctuating stress terms. The 
authors follow Launder23 and model 
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(60) 


1 dp' £ 

- — = 2 c, - a - 
P dx l c k 


c 2 a • V u 

c. 


where c\ c and C 2 c are empirical constants. Comparison with Eq. (39) shows that 
these models have in common the decay of a term and that these would be the 
same if 


c 

K' = 2c, - 
'c k 


(61) 


Besnard, Harlow, and Ra uen zahn^ keep equations for both the turbulent 
heat flux and the quantity A = p'u', which is related to a by 


since they are interested in more complicated equations of state in which the 
relation (19) does not hold. Their a-equation differs from ours in several respects. 
An equation for the Reynolds stress is retained and used in modeling the first 
term on the right-hand side of (40). The triple correlation is broken into two 
terms 


p'u"u" = — 2 p aa + p'u'u' , (63) 

and the latter term is modeled by a gradient diffusion of a. There is a decay of a 
that arises solely from the viscous stress terms. 

IV. SUMMARY AND FUTURE WORK 

We have derived a transport equation for the quantity a, which is the differ- 
ence between the volume- and mass-averaged velocities and is simply related to 
the turbulent heat flux 4 > 6 . Using this equation and an assumption analogous to 
the drift flux approximation of two-phase flow modeling, we have obtained an 
algebraic closure relation for <t>>* that exhibits fluxes due to directed transport 
proportional to — Vp and due to gradient transport proportional to — VT. 

Much work remains to be done before the model can be used in predictive 
calculations of low Mach number flows with large density variations. The equa- 
tion for a involves an additional scalar 6 that is a measure of the density fluctua- 
tions. An equation for b must be derived and terms in it modeled. We hope to use 
the a - and 6 -equations in conjunction with a k-i turbulence model. The k- and in- 
equations must be reexamined to see what modifications are needed when the 
flows have large density variations. When mass transport is important, such as 
in many combustion problems, expressions for the turbulent mass flux must be 
developed. 

In an effort to test some of the modeling assumptions we are currently 
writing a one-dimensional code that solves the turbulence equations of this 
paper. Computed results will be compared with experimental measurements of 
Rayleigh-Taylor instability, turbulent premixed flames, and flows with centri- 
fuging and density variations. These results and extensions of the model will be 
reported in future publications. 
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This research is motivated by the improved ballistic performance of large- 
caliber guns using stick propellant charges. A comprehensive theoretical model 
for predicting the flame-spreading, combustion , and grain deformation phenomena 
of long unslotted stick propellants is presented. The formulation is based upon a 
combined Euler lan-Lagrangian approach to simulate special characteristics of the 
two-phase combustion processes in a cartridge loaded with a bundle of sticks. The 
model considers five separate regions consisting of the internal perforation, the 
solid phase , the external interstitial gas phase, and two limped parameter re- 
gions at either end of the stick bundle. For the external gas-phase region, a set 
of transient one-njimensional fluid-dynamic equations using the Eulerian approach 
is obtained; governing equations for the stick propellants are formulated using 
the Lagrangian approach. The motion of a representative stick is derived by con- 
sidering the forces acting on the entire propellant stick. The instantaneous 
temperature and stress fields in the stick propellant are modeled by considering 
the transient ax i symmetric heat-conduction equation and dynamic structural analy- 
sis. For the internal perforation region, a set of one-dimensional transient 
fluid-dynamic equations is formulated with a coordinate system attached to the 
moving stick. Major distinctions between the present and the conventional formu- 
lations foe interior ballistic simulation are delineated. 


I. Introduction 

Recently, there has been an increasing inter- 
est in the use of stick propellant charges in 
large-caliber gun systems. Stick propellants offer 
many advantages over conventional randomly packed 
multi -perforated granular propellant charges. The 
regular geometry of stick propellants allows a 
higher loading density, flexibility in charge 
design, and easier charge loading. The higher 
charge density is preferable for low vulnerability 
anmunition (LOVA) propellants, which require a 
higher propellant mass to produce an equivalent 
performance. It has been observed by Robbins et 
al . (1,2) that flow resistance through the charge 
of a stick propellant bundle is lower than that 
through packed beds of granular propellants, thus 
enabling faster and more reproducible flame spread- 
ing through the stick propellant charge. The lower 
flow resistance also reduces considerably the phe- 
nomena of high pressure gradients and severe 
pressure waves in the gun tube, as indicated by 
Minor (3) . 

A number of studies (1-10) on various aspects 
of stick propellant combustion have been reported 
to date. The NOVA code, developed by Gough (11,12) 
for ballistic performance of granular charges, was 
used with some modifications (1,4) to predict the 
performance of a stick propellant charge. Results 
obtained were in good agreement with experimental 
data for multiperforated granular NACO propellants 
and single-perforated slotted stick bundles (1). 
However, the same is not true for the case of 
single-perforated, unslotted stick propellants. 

The structure mechanics consideration in the 
continuum modeling of unslotted stick charge com- 
bustion of the modified NOVA code (4) is rather 
etude, due to the application of a steady-state re- 
lationship between radial and hoop stresses and 
internal and external pressures. Although pressure 


distributions in the internal perforation were cal- 
culated, only the external pressure was used in 
evaluating the axial stress component, which is in 
turn related to the intra-granular stress. Grain 
deformation and fracture of unslotted long sticks 
are mainly due to the radial expansion and attain- 
ment of a critical hoop stress. 

Rupture of stick propellants during combustion 
in a gun barrel was observed by Robbins and Horst 
(2) - Grain fracture can lead to high peak pressu- 
res due to increase in total burning surface area. 
The pressure difference across the web of a stick 
propellant can cause the grain to deform prior to 
fracture and alter the flow-channel width and the 
distance between opposite burning surfaces, thereby 
changing the combustion process. Hence, the pred- 
iction of grain deformation and rupture should be 
given due importance in the overall interior bal- 
listic cycle. 

The proposed model is based on a relatively 
new method, previously applied to two-phase react- 
ing flow problems such as spray combustion of 
liquid droplets (14) . In the separated -flow ap- 
proach, the continuous phase (gas phase) is treated 
by a Eulerian approach, while the condensed phase 
media (stick propellants) are divided into several 
representative groups and tracked, using a Lagran- 
gian approach, as they move in the continuous 
phase. The interstitial gas-phase region in the 
present problem is treated in a similar fashion to 
the continuous phase in spray combustion, while the 
internal gas-phase regions in the grain perfora- 
tions are treated separately. The flame- spreading 
and combustion phenomena inside the perforation are 
similar to those in a cylindrical side-burning 
rocket motor grain. The ignition transient analy- 
sis developed by Peretz et al. (15) is therefore 
adapted to model the flame spreading and combustion 
processes in the perforation. This two-phase sep- 
arated flow approach is convenient because the 
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stick propellants have identical dimensions and 
symmetry about their own axes as well as a central 
axis. The number of representative sticks in a 
stick propellant bundle is relatively low, there- 
fore making it possible to study some of the 
detailed flow and combustion phenomena of these 
representative stick propellant grains. 

All of the gun interior ballistic codes and 
analyses developed so far do not consider the dis- 
tribution of chemical species (16-18). The heat 
release is assumed to occur at the same location as 
that of pyrolysis of the propellant. However, for 
LOVA or other types of modern propellants, same 
chemical species could be pyrolyzed from the pro- 
pellant surface, be carried by the flow during the 
ignition transient interval, and then react at some 
downstream location. Therefore, in general, the 
heat release could occur at a different position 
than the site of initial pyrolysis. These complic- 
ated phenomena, which were simplified in previous 
analyses, are included in the present model. 

Objectives of the current research are to for- 
mulate a combined Eulerian-Lagrangian model for the 
combustion of a bundle of stick propellants inside 
a gun chamber. It is intended to cover many as- 
pects of realistic simulation of stick propellant 
combustion characteristics so that the model is ca- 
pable of predicting the phenomena of ignition, 
flame spreading, and combustion of stick propel- 
lants. Development of the model is also intended 
to help explain the experimental observations from 
test firings being conducted by the authors in a 
simulated gun system. The data to be obtained will 
be used for model validation. 

II. Method of Approach 


Physical Model 


Figure 1 shows a bundle of stick propellants 
loaded in the combustion chamber of a gun. For the 
present model, only unslotted long stick propellant 
grains are considered. In the theoretical formul- 
ation, the combustion chamber is divided into five 
separate regions: 1) lurrped parameter region near 
the base pad; 2) internal perforation region; 3) 
external interstitial gas-phase region; 4) solid 
propellant region; and 5) lumped parameter region 
near the base of the projectile. Although the 
solid propellant region consists of many stick pro- 
pellants, only a few representative sticks are 
required to be modeled due to the similarity of the 
sticks in the same family. Each region has a sep- 
arate set of governing equations which are coupled 
through various boundary conditions. The erosive 
burning of the stick propellants under cross-flow 
conditions is also taken into account* Since frac- 
ture phenomena of stick propellants under dynamic 
loading conditions are under investigation, the 
present model is limited to the time period before 
one or more stick propellants rupture as a result 
of high pressure differential across the web of the 
stick propellant. After the onset of rupture, a 
stick propellant could become partially slotted, 
broken, and/or highly deformed. It is then diffi- 
cult to distinguish the internal versus external 
surface. These phenomena are beyond the scope of 
the present model. 


2. Mathematical Formulation 


A. Basic Assumptions: 



Fig.l Stick propellant Charge in a Cartridge 
of a Large Caliber Gun. 


(1) All the stick propellants have the same 
initial geometry and physical conditions. Further- 
more, it is possible to divide the stick bundle 
into a few families so that the calculations can be 
performed for only a few representative sticks. To 
simplify the nethematical formulation, the combus- 
tion of sticks in a bundle is represented by a 
single stick. The rrathematical format can be read- 
ily extended to multiple families of sticks. 

(2) Assumptions used in the gas-phase regions 

are: 

a. no body forces; 

b. bulk viscosity y'is negligible; 

c. Soret and Cufour effects are negligi- 
ble; 

d. gases obey Noble -Abel equation of 

state; 

e. all binary diffusion coefficients are 
equal; 

f. Fick's law of diffusion is valid; 

g. flow is one-dimensional transient 
(properties are uniform in r and 8 di- 
rections) ; 

h. total flow area in the external gas- 
phase region is assumed to be b A 
(Dupuit Forchheimer hypothesis (19) T? 
and 

i. turbulence correlations of flow proper- 
ties in the axial direction are consid- 
ered negligible in comparison with the 
product of their mean flow properties. 
However, turbulence effects in the 
transverse direction are errioeddea in 
the empirical correlations. 

(3) Assumptions used for the solid propellant 

are: 

a. burning in the inner and outer surfaces 
of a stick is axisynmetric; 

b. density of the solid is constant; 

c. any subsurface heat release occurs very 
close to the surface and therefore can 
be lumped onto the surface; 

d. torsion and rotation are negligible; 

e. each stick propellant is locally axi- 
synmetric; 

f. burning at the end surfaces is uniform; 

g. end surfaces are perpendicular to the 
axis of the stick; and 

h. propellant material behaves as a linear 
viscoelastrc material in shear and 
elastic material in bulk deformation. 

B. Overall Structure of the Mathematical 
Model 
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In order to give an overview of the structure 
of cne mathematical model. Fig. 2 was constructed. 
It shows the major components and subcomponents of 
the formulation for each region in the physical 
model. To simplify the diagram, the cross-coupling 
ines between each region are not shown. However, 
it is important to note that the physicochemical 
processes in these regions are closely coupled. 
Mathematical representations of each component and 
subcomponent are given in the following sections. 

C. Governing Equations for the External 
Gas-Phase Region 

The external gas-phase region occupies the 
whole interstitial space in the cartridge, exclud- 
ing the shaded region shown in Fig. 3. By using 
the control volume analysis (reference to Euler i an 
coordinate), the following governing equations are 
obtained . 


CHAMBER WALL 



Fig. 3 Schematic Diagram Showing the Mass 
Fluxes Entering and Leaving the Exter- 
nal Gas-Phase Region (Unshaded and 
Multiply Connected Region) . 



where 7 sin 3 Q is the horizontal component of the 
gastif ication velocity normal to the solid propel- 
lant surface. The second term in the right-hand 
side represents the momentum transfer from solid 
propellant to the external gas-phase region due to 
burning. 


Energy Equation: 

j(o E ) 3 ( o v 1) f ) + > 

% e e e + ge ge e e + ge*e 


3t 

3z 

N 


2 E r 

_ 

j»l °i 

M 

— J — 1 p r u 

n 2 s be 

r y 


’,r ie ie 2 


r (-~^ slnG +*) 

\ k ge ° y 


■ H (T -T ) - 
te ge se 


D u _ p 
ve s e 3t 




' 2 Q /R 

oz w c 


^ere h^ is the enthalpy of the ith species defined 


h i * 4h f,i + j T C p,i dT < 5 > 

o 

^ the sum of the internal and kinetic 
energ? of the gas phase in the control volune. The 
® nGr ^f transfer due to molecular species diffusion 
bulerra” neglecte6 because of a high degree of tur- 

Species Continuity Equation: 

P Y ) 3(<t> p u y ) 

le + e S e 3 e ie * 

* 3z 

3 9 y . 

— (tf> p V ^ Ie ) + (lj > ( 

oz e ge 3z i e v 

"****. “** source term <& , consists of 
contributions cause by surface 1 iyrolysis and gas- 
phase reactions, i.e. ^ 

■ Ou.) + (w ) (7) 

1 e i e,s i e,g K ' 

where the surface pyrolysis part can be expressed 


(tii.) 

i e , s 


r. p A Y * 
be s se is 


* 

where y ig represents the mass fraction of the ith 
species pyrolyzed from the solid propellant before 
mixing with ambient gases. Following the flame 
model proposed by WU et al. (20) in their study of 
erosive burning of homogeneous propellants, the 
solid propellant pyrolyzes into three groups of 
species : 


Solid Propellant 


Oxidizer-rich gases (0) 


Fuel-rich gases (F) 

+- 

First group of species 
with delayed reaction fhR' 
( 9) 

Under low cross flow conditions, the flame struc - 
ture adjacent to a burning homogeneous propellant 
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surface exposed to a large cavity can be assumed to 
nave three stages, as shown m Fig. 4a. This flame 
structure is based upon the erosive burning study 
in a rocket motor by Wu et al. (20). In the case 
of large-caliber guns densely loaded with stick 
propellants, the void spaces adjacent to burning 
surfaces are relatively small , and the species py - 
rolyzed from the surface can be entrained by the 
high-velocity gases flowing along the axis. The 
heat release in the final flame generated by the 
chemical reaction of pyrolyzed species from a spe- 
cific location occurs at a downstream location, as 
shown in Fig. 4b. To determine the gas phase reac- 
tion rate and heat-release rate, the same chemical 
reaction mechanism, proposed in Ref. 20, is adop- 
ted. This mechanism can be represented by three 
overall chemical steps: 
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F + v 0 
F 0 


DR1 
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DR2 
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Oxidizer-rich gases (0) can be regarded as 
species, fuel-rich gases (F) as a group of aldehy- 
des (CH~0) , and other species generated by surface 
pyrolysis as one group delayed reaction species 
(DRl) . Based upon the work of Fifer (21) and 
Kubota (22), chemical reaction in the Fizz zone is 
largely due to the reaction involving reduction of 
NO to NO. After the delay in the dark zone, the 
reactions in the final flame can be assumed to 
occur at- high activation energies, associated with 
reactions (11) and (12) to form the final products 
from DRl and the second group of delayed reaction 
species (DR2) . In the dark and final flame zones, 
chemical reactions result in oxidation of CO, and 
perhaps some by NO. More detailed discussions 
of kinetic parameters and mechanisms are given in 
Ref. 20-22. 


It is useful to point out that the chemical 
reaction mechanism proposed in Ref. 20 can be fol- 
lowed for highly turbulent cross-flow situations. 
The 0 and F species pyrolyzed from the solid pro- 
pellant originate at the same place, and flow 
together in a torturous path to form DR2 and final 
products (P) . The delayed reaction species of DRl 
can also be considered to flow together in the 
process to form product species P. Gases in the 
cross-flow can transfer heat to these species and 
alter the delay times required to form product spe- 
cies P. The stoichiometric coefficients V v n , 
' and v are therefore functions of prope r - 
lane ingredients only. The method for calculating 
these parameters is based upon the original molecu- 
lar structure of the propellant and is discussed in 
Ref. 20. 


In Ref. 20, the rate of production of species 
"i" was based upon the chemical reaction cate as 
well as the eddy-break-up rate (23) controlled by 
the turbulence intensity and concentration gra- 
dient. Turbulence intensity is so high in the gun 
situation because of base-pad ignition and combus- 
tion of propellants that the species diffusion term 
can be regarded as extremely short. Therefore, the 
rate of consumption or production of species is de- 
termined solely from chemical reaction time. 


The rate of chemical reaction of species F, 
DRl, 0, and DR2 are given below in the same form as 
m Ref. 20. 
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Fig. 4 Flame Structure of a Homogeneous Solid 
Propellant (a) Two-dimensional Struc- 
ture under Low Cross-Flow 

Velocities, (b) Distorted Flame under 
Extermely High Cross-Flow Velocities 
(Nearly One-Dimensional Structure). 


Equation of State: 

The Noble-Abel dense gas law is used. 

p (— b) = RT (17 > 

e ge 
§e 

The initial and boundary conditions as well as em- 
pirical correlations for the external gas-phase 
region are given in a later section. 

D. Governing Equations for the Stick 
Propellant 

The following governing equations for a repre- 
sentative stick are derived based upon the 
Lagrangian coordinate. 

1) Mass and Momentum Equations: 

The instantaneous mass of the representative 
stick can be calculated from 
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M (t) 
s 


[L(t)/2 

D TT[r 2 (t,D - r 2 (t,5)ld? (18) 

J-L(t)/2 S ° 1 


The instantaneous values of stick propellant 
length, L(t) , and the local inner and outer radii, 
r (t,S) and r (t,£) , are determined by integrating 
tAe following°first order differential equations. 


3r i (t.C) 

3t 

3r (t,0 
o 

3 1 


+ V (t,£) 
bi sr^ 


be 


(t + v gr (t,Q 


(19) 

( 20 ) 


2) Transient Heat Conduction Equation: 


To determine the instantaneous temperature 
distribution in the stick propellant, a heat con- 
duction equation must be considered. The equation, 
which takes into account the subsurface radiation 
absorption for translucent propellants, has the 
following form. 


small 


3(C T ) . . 9T 

9 s . I ±_ ( r k — §.) + 
Dt r 3r U s 3r J 



■["* 


-v * < w 


(24) 


^■-' b (t > + V SRB (t) - V SLB (t) (21 > 

where V and V represent the radial velocities 
of the sr Inner a&S°outer surfaces of the stick pro- 
pellant due to mechanical deformation with respect 
to the centerline of the stick. and re- 

present the rate of mechanical deformation or the 
right and left boundary surfaces with respect to 
the geometric center of the stick. The equation of 
motion is formulated according to the following 
momentum balance principle. 

d(M u ) 

— — — « E + (net rate of momentum flux flowing 
into the control volume encompass- 
ing the stick propellant) 

( 22 ) 

which gives (see Fig. 5; 


where E, i| the black-body emissive power evaluated 
as oT , and "a" represents the flux absorption 
coefficient of the propellant. This equation is 
based upon a two-flux model which assumes that the 
radiation fluxes are dominant in the radial direc- 
tions (inward and outward) . The source terms 
represent the net rate of energy absorbed due to 
radiant energy fluxes. The outward and inward 
radiant fluxes, I and J , can be determined from 
the following flux-transport equations. 

d(rl ) 1 

— = - ( s+a ) r I r + arE b + + - sr(I r +J r ) (25) 

d(rJ ) 1 

dr - = (s+a)rJ r - arF, b + J r " y sr ^ I r +J r ^ 

These two equations were used by Gosman and Lock- 
wood . ( 24 ) 
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where F represents the net force acting between 
the stick propellant and chamber wall, and F is 
the force between adjacent propellants. ^ 


3) Dynamic Structural Analysis 

As a result of the different ignition and 
flame spreading processes occur ing in the internal 
perforation and exterral gas-phase regions, a pres- 
sure differential exists across the web of a 
representative stick. A finite element analysis is 
needed to compute the resulting transient viscoe- 
lastic deformation of the stick propellant and to 
predict the attainment of a critical condition for 
grain fracture. Regression of the boundary as a 
result of pyrolysis and burning should also be 
taken into account. Furthermore, the mechanical 
properties of the stick propellant must be specif- 
ied. 


The propellant material can be considered as 
linear viscoelastic in shear and elastic in bulk 
deformation. This is a conmonly accepted practice 
for solid propellant (25, 26). The elastic bulk 
behavior is assumed to follow 


a 


i 

i 



(27) 



where K is the bulk modulus. The deviatonc beha- 
vior is taken as 


3e ( t ’) 

dt’ (28) 


where the shear relaxation modulus G,(t) is assumed 
to be of the form 

Fig. 5 Momentum fluxes and Pressure Forces 

Acting on the Entire Surface of a G (t) = G + (G -G (29) 

Single Perforated Stick Propellant. 1 00 ° “ 
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(33J 


where is the long-time shear modulus, G is the 
short-time shear modulus, and S is the decay con- 
stant. Since a closed- form solution for the 
axisynmetric dynamic problem is not possible, a 
well-established finite element code "HONDO- II" 
(26) is employed for computations of grain deform- 
ation. The code utilizes the principle of virtual 
work for the solution. It states that at all the 
points along the path of motion, the differential 
virtual work, 6tt r must vanish for all variations 6x 
satisfying the imposed displacement boundary condi- 
tions (26) . <5 tt i s defined as 


6 TT 


Px 5x^dv + 


O^Sx, _dv - | pf^iSx^dv 


k ,ra 


- 0 T <5x_ ds 
k 


(30) 


where T*is surface traction and a ^ is Cauchy 
stress tensor. The stick propellant is divided 
into a number of elements over the cross-sectional 
area of the web. The HONDO- I I code uses four mode 
bilinear isoparametric elements. The basic equ- 
ation of motion, viz., the minimization of virtual 
work, is then considered for each of the nodes in 
the stick propellant. Thus, the equation of motion 
for a node becomes 



"k a* „ 
px 9, dv + 
~k 


km , 

0 6, dv 

,ra 




pf 6, dv - (pi S^ds 


k„St, 


- 0 


where 9^ = [$*♦ $ 4>£] 


(31) 


l 

and are the bilinear interpolation functions. 
In Eq. (31), W is the total number of elements sur- 
rounding the node in question. The calculations 
are carried out on an element -by-element process to 
get the final equations of motion. The time in- 
tegration of these equations gives the positions of 
nodes at the new time step. A central difference 
method is used for time integration in the HONDO-II 
code. 


E. Governing Equations for the Internal 
Gas -Phase Region 

Conservation equations for the control volume 
in the perforation of a stick propellant are simi- 
lar to those given in Ref. 15. 


Continuity Equation: 
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gl gi 9 

Momentum Conservation Equation: 
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Energy Conservation Equation: 


The energy conservation equation written for 
the stored total energy (internal and kinetic) per 
unit mass, E . , is 

gi 
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Species Continuity Equation: 


The species continuity equation for the inter- 
nal gas-phase is similar to that for the external 
gas-phase, with the void fraction taken as one and 
the absolute velocity replaced by the relative ve- 
locity. 


3[p 

di + 


35 



+ ( Vi 


(36) 


small 


where the source term (w. )., similar to ( ) in 
Eq. (7) , consists of contributions due to surface 
pyrolysis and gas-phase reactions. The flame model 
for the internal gas-phase region is the same as 
that given earlier for the external gas-phase 
region. 


F. Heat Losses to the Walls of the Combustion 
Chamber 


In order to consider heat losses from the com- 
bustion zone to the cartridge chamber, the gun 
tube, and the projectile, temperature profiles in 
these metal components are required. For a test rig 
with a blowout diphragm and short barrel, the tran- 
sient ignition and combustion phenomena occur in an 
extremely short time interval and the heat loss to 
the surrounding walls can be considered negligible. 
However, if the simulation is made for the full 
ballistic cycle occurring in a long gun barrel, the 
loss to the walls must be considered. The heat con- 
duction to the gun tube can be given as 


3(C T ) 
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r 3r w 


ai 
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(37) 


where £ is the distance measured from the center of 
the stick, and V . is the velocity of the gas rela- 
tive to the velocity of the stick propellant. V ■ 
is related to the absolute velocity u ^ by 


The transient heat conduction equation for the 
wall on the breech end will follow the one- 
dimensional form: 
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3(C T ) ~ 3T 

w vr 3 , , w N 

p = — ( k ) 

K w dt dz w dz 


2) • Boundary Conditions for the External and 
(38) Internal Gas-Phase Regions 


The equation for the projectile is: 


D(C T ) . 3T 

w w 3 /t vi 

P w Dt dz w 3z 


(39) 


The solution of these equations are coupled to 
the external gas-phase region through the boundary 
conditions on the wall surfaces. 


For simplicity in expressing the boundary 
treatments, only one representative group of stick 
propellants is considered, then each stick propel- 
lant has the identical instantaneous velocity and 
length. Therefore, the governing equations for the 
left lumped-parameter region (see Fig. 6) or called 
left control volume (LCV) can be readily derived 
and written as follows. 


G. Initial and Boundary Conditions 

For each of the above governig equations, 
there is a set of boundary and/or initial condi- 
tions required to complete the formulation. The 
initial conditions can be specified readily, since 
the gas and stick propellant velocities are zero 
and the pressure and temperature of gas in the car- 
tridge is at room conditions. The stick propellant 
charge with known geometry and its surrounding 
chamber are also at room temperature. Since the 
amount of .air in the initial loading of the car- 
tridge is extremely small in comparison with the 
gases generated from combustion, it can be treated 
as any one of the five species discussed above, in 
view of the fact that air contains oxygen, it is 
treated as oxidizer-rich gas (0) . 


A number of important boundary conditions are 
given in the following. 

1) . Boundary Conditions for the stick Pro- 
pellant 


The boundary condition for Eq.(24) at outer 
surface of the stick propellant can be written as 
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where ^PIP J e re P resents the rate of energy input 
due the deposit of condensed ph$se igniter products 
onto the external surface and (q .) the net ra- 
diative heat flux to the surfaci cln be expressed 
as 
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Similarly, the boundary condition for the inner 
surface of the stick propellant can be written as 
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The pressure distributions along the internal per- 
foration and external interstitial regions solved 
from the gas-phase equations are used as boundary 
conditions • for the solid-phase dynamic structral 
analysis. 


Mass balance: 
dp^ 


g 


LCV 


dt 


= - P„ 


(u +r ) 
S b LB 




'LCV L 

P o 

8 lcv 


s p . LCV 
ign 


m -m -m 
lgn oute out! 


(44) 


7TR <J> Z 
c LCV L 


where z^ is the length of the left lumped -parameter 
region and can be calculated from 
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Momentum balance: 
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Fig. 6 Control Volume Considered for Deriving 
Governing Equations in the Left 
Lumped -Parameter Region. 


CRJGSNAL PAGE 15 

OF. POOR QUALITY 



Energy balance: 


region remains the same as that in the first case. 
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Equation of state: 
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The void fraction in the left lunped -parameter 
region, T can be evaluated from the continuity 
equation or the condensed-phase products of the lg- 
niter which is given as follows. 
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In the third case, the flow velocities at the 
left boundray of both the perforation region and 
the interstitial regions are negative. Therefore, 
two additional compatibility relationships are 
available. In this case, both the internal and ex- 
ternal gas-phase boundaries are treated in a 
similar fashion as that for the internal perfor- 
ation region in the second case. 

The void region between the right end of the 
stick propellants and the projectile base is a 
lumped -parameter region, which could contain either 
a base pad with powders of igniter material or an 
ullage space. In the case of base pad, a set of 
equations similar to Eqs. (44) - (51) are derived. The 
right control volume (RCV) is considered to move 
with the projectile and the left surface of RCV 
coincides with the right boundary of the stick 
bundle. A detailed discussion on the compatibility 
relation mentioned above is given in Ref. 27-28. 

H. Empirical Correlations and Constants 

Several empirical correlations are used in the 
model. For the internal gas-phase region, these 
comprise of the correlations for the burning rate 
law including erosive burning effect, the convec- 
tive heat-transfer coefficients, and viscous drag 
coefficient. The latter two correlations are same 
as those used in Ref. 15. New erosive burning cor- 
relation being developed by authors and coworkers 
in parallel to this study will be used. For the ex- 
ternal gas-phase region, correlations for the 
convective heat-transfer coefficient, and viscous 
drag coefficient between the gas-phase and solid 
surfaces are needed. These correlations are similar 
to those for the internal perforation region except 
the fact that they are based on a hydraulic diame- 
ter defined as 


(51) 

Ttie boundary values of velocity, density, 
pressure, and temperature on the left hand side of 
both internal and external gas-phase regions can be 
solved from a number of relationships coupled with 
Eqs.(44) and (48) -(51). Depending upon the flow di- 
rections at left boundary of the internal and 
external regions, three different cases are identi- 
fied. In the first case, the flow velocities at 
stick perforation and the interstitial regions are 
both in the positive z direction, it can be assumed 
that the pressure at LCV is equal to those at the 
left boundary of the stick perforation and the ex- 
ternal gas-phase region. Also, the stagnation 
temperature of the gas in LCV can be assumed to be 
the same as those on the left boundary. One com- 
patibility relation along the characteristic line 
in each gas-phase region is used together with the 
aforementioned equations for solving the flow prop- 
erties in the LCV and the left boundary. 


2( R 2 - Nr 2 ) 
n = c o 

H R + Nr (52) 
c o 

The burning rate law for the external region 
will be the same as that for the internal perfor- 
ation region. 

For the stick propellants, correlations for 
F pn and F pw' the criterion, and mechanical 
ana therms! properties are needed. At the present 
time, not all of this information is available; es- 
pecially F^ and need to be characterized. 

4. Summary of Differe nces Between the 
Present and Conventional Formulation 

In order to bring out the major differences 
between the present formulation and the conventio- 
nal interior ballistic predictive models (4), a 
suirmary table is given below. 


In the second case, the flow velocity at the 
left opening if the stick perforation is negative, 
but the bulk velocity at left boundary of the in- 
terstitial void region is still positive. In this 
case, one additional compatibility relationship can 
be used .for the stick perforation region. Only the 
pressure at the opening of the stick perforation is 
considered to be equal to that in the LCV. The 
treatment for the left boundary of the interstitial 


As suirmerized by Robbins and Einstein (29), 
there are differences between the measured and cal- 
culated pressures from the NOVA code or its 
extensions for long unslotted stick propellants. 
Also, the measured muzzle velocities are higher 
than those calculated for slotted stick propel- 
lants. A number of improvements and considerations 
suggested in the workshop (29) are incorporated in 
the present model. 
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Differences Between the Present Formulation ana the 
Conventional Interior Ballistic Formulation 


Subject under 
Consideration 

i 

Present Formulation 

Convential formulation (4) 

Typical grain 
configuration 

* Simulation of a number of 
typical full-length grains 
in a bundle of stick prop- 
ellants. 

* Each stick is modeled as 
separate tube with deform- 
able and combustible wails. 

* Simulation of an average 
grain in a spatial location 
along a packed bed of gran- 
ular propellants. 

* Each bundle is modeled as a } 

continuum characterized by j 

the velocity and stress in j 

the sticks . | 

Averaging of 
flow 

properties 

* The external flow propert- 
ies are averaged over the 
cross-sectional flow area 
of interstitial voids, 
while the internal flow 
properties are averaged 
over the flow area of each 
stick propellant. 

* Flow properties in both ext- 
ernal and internal regions j 

are averaged over their res- 
pective flow areas. 

i 

j 

Grain 

deformation 
and fracture 

* Simulated by the unbalanced 
pressure forces between the 
internal perforation and 
external interstitial void 
region. 

* Linear viscoelastic consti- 
tutive law is used. 

* Employs dynamic finite- 
element structure mechanics 
computational code. 

* The process of grain defo- 
rmation and fracture are 
not addressed, except the 
longitudinal stresses are 
considered in the solid- 
phase momentum equation. 

* Linear elastic constitutive 
law is used. 

* Employs a steady-state rel- 
ationship between stresses 
(radial and hoop) and press- 
ures (internal and external) j 

! Grain 

displacement 

and 

acceleration 

* The kinematics of the full- 
length grain is determined 
from the summation of all 
forces exerted on the grain, 

* The bulk properties of the 
grains are determined from 
local momentum balance. 

j 

Radiative heat 
transfer 

* Subsurface radiation pene- 
\ tration is allowed and tre- 
' ated by a two-flux model. 

1 

* No subsurface radiation ] 

penetration is considered. 

: Type of 
• formulation 
and frame of 
reference 

i 

* Kinematics and gram defor- 
mation are formulated by 
following the stick (Lagra- 
ngian approach) while the 
gas-phase properties for 
internal and external regi- 
ons are determined from a 
fixed frame of reference 
(Euler lan approach) . 

* Both the gas-phase and solid 
-phase properties are all 
determined from the conserv- 
ation equations formulated 
based upon a fixed frame of 
reference (Eulerian 
approach) . 

Species 
distribution 
and location o 
heat release 

* Five groups of species are 
considered. 

t 

* Heat release does not have 
to occur at the site of 
pyrolysis. 

* Gas phase is made of combus- 
tion products from ignition 
and propellant. 

* Heat release occurs locally 
at the site of pyrolysis. 
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Letter Symbols 


Nomenclature 


a 

A 

A. 

l 


flux model absorption coefficient, nf 1 
cross-sectional area of the gun barrel, m 2 
preexponential factor of the 1 species 
m /kmol-s 

cross-sectional area of the perforation, m 2 


Ag specific area of the external surfaces of 

the stick propellant, m 1 

b covolume of the Noble-Abel equation of 

state, m Ag 

C p constant pressure specific heat, JAg-K 

C s specific heat of stick propellant, JAg-K 


constant volume specific heat, JAg-K 


DR1 group of species pyrolyzed from propellant 
surface having delayed reactions 
DR2 delayed reaction species generated from 0 
and F species 

D y viscous drag force per unit area, N/m 2 

' binary diffusion coefficient, m 2 /s 

e lj deviatonc strain 

E total stored energy (internal plus kinetic) 

per unit mass, JAg 
E a activation energy, J/kmol 


Ej^ black-body emissive power, = ;T 4 , J/m 2 -s 

f body force per unit volume in k*'* 1 direction. 

N/nT 

' fuel rich species pyrolyzed from propellant 

external force exerted on the solid propel- 
lant in the axial direction, N 
h specific enthalpy, JAg 

h c convective heat-transfer coefficient, W/in 2 -K 

hj specific enthalpy of the j th species, JAg 

h t total heat-transfer coefficient, W/m 2 -K 

o 

Ahf ^ standard enthalpy of formation of the i fc ^ 

' species, JAg 


outward radiation flux in the positive 
radial direction, W/m 2 

inward radiation flux in the negative radial 

direction, W/m 2 

thermal conductivity, W/m-K 

bulk modulus of the propellant material, 

N/nT 

total mass of gas in control volmve, kg 


0 

P 

P 

P 

b 


instantaneous total mass of a single stick 
propellant, kg 

oxidizer species pyrolyzed from propellant 

pressure , N/m 

final product species 

perimeter of the internal perforation, m 


q rad radiativ e heat fulx per unit time absorbed 
by the solid propellant surface, W/m 2 
Q s surface heat release due to pyrolysis, JAg 
chem 


rate of heat loss to the tube wall, W/m 2 


r b propellant burning rate, m/s 

r^ inner radius of the perforation, m 

r Q outer radius of the stick propellant, m 

R gas constant, JAg-K 

R c radius of the combustion chamber, m 

R y universal gas constant, J/kmol-K 

s flux-model scattering coefficient, 1/m 

deviatoric stress tensor, N/m 2 

5 surface, m 2 

t time, s 

T temperature, K 

Tf adiabatic flame temperature of the solid 

propellant, K 

u absolute velocity, m/s 

V gas velocity relative to the solid propel - 
y lant, m/s. 

^ volume, rn 

molecular weight of the i th species, kg/kmol 
\ coordinate axis in k tn direction, m 

x* k acceleration in k th direction, m/s 2 

^ X k,m deformation gradient tensor 

Y i mass fraction of i th species, i could 

represent F, 0, DR1 , DR2, or P 
z axial coordinate in Eulerian coordinate 
system, m 


kk 


0 . 


km 


w 

U)< 


Greek Symbols 
dilatory strain 

surface emissivity of solid propellant 

angle measured ocw from axis to the 
tangent to the inner surface of the stick 
propellant, rad 

angle measured ccw from axis to the tan- 
gent to the outer surface of the stick 
propellant, rad 

dynamic viscosity of the gas, N-s/m 2 
number of kmoles of i species 

Lagrangian axial coordinate, m 
virtual work, J 
density, kg/nr 

Stef an-Bol tzmann constant, w/m 2 -K 4 
stress tensor, N/m 2 

viscous shear stress, N/m 2 

surface traction in k th direction, N/m 2 

tube wall shear stress, N/m 2 

normal viscous stress. N/m 2 

mass fraction of ccmbustion product in the 

condensed phase from the igniter 

rate of production of i species, kg/m 3 -s 


Subscripts 

e external interstitial region 

g gas-phase, internal or external region 

i internal perforation region 
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_ :n igniter 

left ooundary of the stick propellant bardie 
LCV left lumped -parameter region in the car- 
tridge 

RB right boundary of the stick propellant 
bundle 

s solid propellant or surface 
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Editor's Narrative: 


Round-Table Discussion: Propulsion Applications of Mixing and 
Demixing Processes of Multiphase Flows 

Several points were made either explicitly or implicitly in 
the presentations and discussions on the second day concerning 
mixing processes in multiphase systems and the modeling of these 
systems. First, there are situations where transport in turbu- 
lent flows can be in a direction other than along the gradient of 
the property being transported. This comes about through cou- 
pling of density inhomogenieties with pressure gradients, the 
former being a consequence of the large temperature differences 
in combustion environments. Hence, in turbulent flames, some 
transport processes are "counter gradient" or more correctly, 
"non-gradient" in nature. Errors due to neglect of these pro- 
cesses are not understood, but could possibly be significant in 
some combustion systems. Modelers are in the process of learning 
to incorporate these ideas into practical codes. 

Second, the weakest element in modeling of turbulent multi- 
phase flows is probably the turbulence modeling itself. The 
treatment of the coupling between the turbulence and the non- 
continuously distributed dispersed phase is also a critical area 
and includes the question of whether, for numerical purposes, to 
t- re 3t the dispersed particulate (or droplet) phase as a continuum 
or as an assembly of separate particles. 


PRKCECP 1 ’" 1 


ANX 


f YjS F&MED 


195 



Third, it is possible to treat a very complex problem in 
great detail by using appropriate levels of approximation and by 
coupling separate regions together which are best treated with 
different models or with different levels of approximation. 

Fourth, while continuum approaches to multiphase flow calcu- 
lations are not particularly illuminating in revealing details of 
the physical processes at work (due to the necessity of averaging 
quantities before solving the governing equations) , there are 
cases where they can provide useful results for a given problem. 
The two-fluid approach works best where the dispersed particulate 
phase is monodisperse, and where the flows are non-reacting. 

Sixth, advantages and disadvantages of Favre averaging 
approaches were discussed. An advantage is that the conservation 
equations for variable density with Favre averaging are much like 
the standard constant density equations. Disadvantages are that 
there is some difficulty in obtaining molecular terms and com- 
puted and measured quantities are more difficult to compare. 

Finally, physical situations for which multiphase flow mod- 
els are appropriate generally contain a very broad range of 
length and time scales within the same problem. This means that 
when computational approaches to these problems are used, some 
levels of approximation will need to be made in order to deal 
with these large ranges of scale. The point in the analysis at 
which approximation is made is a key difference between a two- 
fluid model and an Eulerian-Lagrangian model of a dispersed par- 
ticulate phase in a fluid phase. 
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